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We investigate a number of fermionic condensate phases on the honeycomb lattice, to determine whether 
topological defects (vortices and edges) in these phases can support bound states with zero energy. We argue that 
topological zero modes bound to vortices and at edges are not only connected, but should in fact be identified. 
Recently, it has been shown that the simplest s-wave superconducting state for the Dirac fermion approximation 
of the honeycomb lattice at precisely half filling, supports zero modes inside the cores of vortices (P. Ghaemi and 
F. Wilczek, 2007). We find that within the continuum Dirac theory the zero modes are not unique neither to this 
phase, nor to half filling. In addition, we find the exact wavefunctions for vortex bound zero modes, as well as 
the complete edge state spectrum of the phases we discuss. The zero modes in all the phases we examine have 
even-numbered degeneracy, and as such pairs of any Majorana modes are simply equivalent to one ordinary 
fermion. As a result, contrary to bound state zero modes in p x + ip y superconductors, vortices here do not 
exhibit non-Abelian exchange statistics. The zero modes in the pure Dirac theory are seemingly topologically 
protected by the effective low energy symmetry of the theory, yet on the original honeycomb lattice model these 
zero modes are split, by explicit breaking of the effective low energy symmetry. 

PACS numbers: 74.20.Rp, 05.30.Pr, 03.67.Lx 



I. INTRODUCTION 

In recent years, p x + ip y fermionic condensate states, 
have received much attention due to the expectation that 
vortices in this state will exhibit non-Abelian (Braiding) 
statistic ;) 1 ! 2 ! 3 ! 4 ! 5 ! 6 ! and their potential applicability to topolog- 
ical quantum computing 78 . A frantic experimental effort to 
observe these effects is currently under way^H The non- 
Abelian effects are caused by the presence of quasiparticle 
zero modes (states w ith energy pre cisely at the Fermi level) 
bound to vortex core d 1 ! 2 ! 3 ! 4 ! 5 ! 1 i | i2 | i3 [T4] Th ese zero modes ap- 
pear at sample edges as wel l 1 ! 15 ! 16 -!, and we will refer to them 
collectively as topological zero modes. 

Zero modes in any BCS mean field Hamiltonian 1 2 3 can al- 
ways be expressed as Majorana fermions. Pairs of Majorana 
states will combine to form single fermionic degrees of free- 
dom, which can then be occupied or not. The p x + ip y super- 
conducting states allow single zero modes bound to vortices 
(of unit vorticity). Therefore, fermionic states can only be 
formed by a superposition of two zero modes, bound to dif- 
ferent vortices. In this way, these fermionic modes provide 
a natural entangled state between two spatially separated ob- 
jects (the vortices) 2 . This entanglement is the source of the 
non-Abelian mutual statistics, when moving one vortex adia- 
batically around another. 

The apparent rarity of p x + ip y superconducting states has 
made difficult the effort to observe these zero modes in exper- 
iment. It would therefore prove useful to have further candi- 
date states for displaying non-Abelian statistics, which could 
then be searched for experimentally. 

The topological zero modes in the p x + ip y state, are found 
as solutions of a set of coupled Dirac-like Bogoliubov-de- 
Gennes (BdG) equations. The source of the Dirac-like be- 
havior is the symmetry of the p x + ip y superconducting order 
parameter. An alternative way to end up with BdG equations 
of the form of a Dirac equation, one which does not require the 
pairing function to be of the p x + ip y form, is to have a kinetic 



energy term that is of the Dirac form. The most celebrated 
example where this occurs is in the honeycomb lattice tight- 
binding model, where close to half filling the band structure 
has Dirac-like dispersion, and the behavior of the system can 
be approximated by two flavors of Dirac fermions. The effec- 
tive Dirac like dispersion has been experimentally observed 
in monolayer graphene 17 . With the Dirac-like behavior of the 
BdG equations already guaranteed in this approximation, we 
are now free to ask whether zero modes exist in vortex cores of 
a whole variety of superconducting states on the honeycomb 
lattice. The simplest state one could consider is the s-wave 
spin-singlet pairing state. Some time ago 18 4 it was shown that 
this same superconducting state in a (square) lattice model 
with Dirac dispersion near the Fermi energy supported zero 
modes bound to vortex cores. More recently 19 zero modes 
were shown to exist (bound to vortex cores) in this state, in 
the Dirac continuum theory of the honeycomb lattice at pre- 
cisely half filling. These findings are in stark contrast to the 
behavior of two dimensional s-wave superconductor vortices 
in fermionic systems with simple quadratic dispersion, where 
no zero modes exist^Sl. 

Following this radically different result, in this article, we 
will investigate other geometries and phases for possible pres- 
ence of topological zero modes. In addition, we will de- 
termine whether the zero modes appear in the actual lattice 
model. 

Even the simplest effective attraction between fermions can 
cause a superconducting state to appear. However, the pre- 
cise nature of the phase, namely the symmetry of the order 
parameter, depends on the details of the effective interaction. 
In Ref. 21 it was shown in a simple mean field analysis that 
fermions on the honeycomb lattice paired in spin singlets may 
support not only an s-wave state, but also an effective p x + ip y 
state, as well as a mixed s-wave/p^ + ip y phase (and earlier 
workp21 a i so suggested a p-wave superconducting state may 
appear in graphene). Given the evidence of zero modes in 
the s-wave phase at half filling^, it is interesting to explore 
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whether the p x + ip y spin-singlet phase may also support zero 
modes, as well as whether these zero modes are peculiar to the 
half filling point (recently^ it was shown that s-wave super- 
conductivity is far more likely to appear in a fermionic sys- 
tem in the honeycomb lattice away from half filling). We will 
show in this manuscript that zero modes appear in the Dirac 
continuum theory in both the s-wave and the p x + ip y phase, 
even when deviating from the special half filling point (The 
Px + ipy state is in fact gapped only when the fermions are 
away from half filling). We also demonstrate that the zero 
modes we find are 4-fold degenerate, and so fermionic modes 
can be formed by pairs on the same vortex. The mechanism 
for entanglement between vortices is therefore unfortunately 
lost, and no non-Abelian effects are expected in these sys- 
tems. Note that it is not obvious that no condensate phase of 
fermions on the honeycomb lattice will exhibit non-Abelian 
statistics. 

As mentioned above, in the p x + ip y state (with the regu- 
lar quadratic kinetic energy) it is known that the zero modes 
bound to vortex cores appear in conjunction with edge states 
zero modes. We argue in this manuscript that this connection 
is in fact quite general, and that vortex core topological zero 
modes and edge state zero modes should in fact be identified. 
We demonstrate this in the superconducting states we analyze 
here, by finding the low energy edge states of each phase. As 
expected, we find precise correspondence with the vortex core 
bound states - in both spin-singlet phases, there exist 4 zero 
modes. A further signature of the identification of the vortex 
and edge states is that the wavefunctions have all the same 
physical parameters - the same exponential decay lengths, as 
well the same oscillation length scales (when those exist). 

A perhaps simpler indication of whether zero modes can ap- 
pear at the edges or vortex cores of a given superconducting 
(SC) state, is to consider the SNS (superconducting-normal- 
superconducting) junction 24 , with some phase difference <fi 
between the SC droplets. We find that as in the regular p x +ip y 
stated zero modes appear only when (p = n. As we will 
see, the edge state calculation we employ in the continuum 
limit is limited in the type of honeycomb lattice edge it can 
be used for. For this reason the SNS junction calculation is 
useful - it shows that at least within the continuum limit the 
precise alignment of the edge is immaterial. We will find that 
the number of zero modes found in the SNS geometry is 8 
rather than 4, giving us the first hint that details of boundary 
conditions are important in this problem - the SNS junction 
geometry has extra symmetries, compared with the edge and 
vortex cases. 

The zero modes we will uncover in what follows, all have 
an even degeneracy. In general (for unit vorticity), only single 
SC quasiparticle zero modes are topologically protected to all 
possible perturbations, however, if there is a symmetry man- 
dated degeneracy, which is not broken by any of the pertur- 
bations, a degenerate set of zero modes can still be protected 
to perturbations, and hence topologically protected (modulo 
symmetry mandated degeneracy). We begin the main body of 
our manuscript with a discussion of this distinction in general 
settings. 

In our case, the 4-fold degeneracy is mandated by the sym- 



metries of the Dirac continuum theory. However, the full sym- 
metry is not present in the underlying lattice model, and so 
there is a danger that the zero modes can appear in the contin- 
uum model, but not in the lattice model. By numerical diago- 
nalization of the lattice model, we find indeed this is the case 
- the zero modes do not exist in the lattice model. 

The remainder of our manuscript is organized as follows. 
In section |n] we discuss topological protection of zero modes 
in conjunction with symmetry mandated degeneracy. In sec- 
tion |lll] we present our general argument for identifying zero 
modes bound to vortex cores and at sample edges. We then 



proceed to section IV where we present the honeycomb lat- 



tice BCS model, describe a number of possible superconduct- 
ing states, and then set up a continuum limit for the model, 
which allows us to perform explicit calculations looking for 
zero modes in the SNS junction (section [7), the edge states 
(section |Vl| i, and finally in the vortex cores (section VII i. We 
present the results of a numerical calculation on the precise 



honeycomb lattice BCS model in section VIII In section IX 



we proceed to discuss a possible experimental realization of 
superconducting states on the honeycomb lattice. We then dis- 
cuss Zeeman splitting in section[Xj and propose how it can be 
used to test the physics we describe here experimentally, us- 
ing the absorption spectrum of the system. We conclude our 
manuscript with the discussion in section [XI] 



II. ZERO MODES IN BCS HAMILTONIANS MODULO 
SYMMETRY MANDATED DEGENERACIES 

The most general BCS 25 Hamiltonian has the form 



ab 



(1) 



where a, b are generalized coordinates, and may include spin, 
position and any other degree of freedom one can imagine. 
The fermionic operators f a satisfy standard anticommuta- 
tion relations. The operator h is hermitian, and the pairing 
function must be anti-symmetric in the generalized coordi- 
nates A ao = —Ah a . A Bogoliubov transformation 7^ = 
J2 a [ u Eafl +VEaf a ] diagonalizes this quadratic Hamilto- 
nian. From the eigenstate equations 



E~/l the 



Bogoliubov-de-Gennes (BdG) equations are derived 



(2) 



where ip = (u a ,v a ) (we drop the E index to avoid clut- 
ter). The BdG equations are then simply the eigenvalue prob- 
lem 7iBdG^ = Eip, in terms of the BdG Hamiltonian 7Ybc1G- 
This operator has the symmetry lo x, H^qijj x = — H^ dG , where 
uj x is the x-Pauli matrix in the so-called Nambu spinor basis 
w I,!, ' z (see table acting on the (u, v) components of HedG- 
The Nambu spinor obeying ui z ip ~ +ip corresponds to a pure 
fermion (v a = 0), while ui z ip — —ip corresponds to a pure 
hole (u a — 0). This relation tells us that given an eigenvec- 
tor ip with eigenvalue E, u> x ip* will also be an eigenvector 
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a,p 


(1,V 


A, B 


index values 




U 


1,2 


R, L 



TABLE I: Pauli matrix sets 



with the eigenvalue — E. The matrix eigenvalue equation (|2]i 
can indeed yield both positive and negative eigenvalues, how- 
ever it is important to note that in this formulation of the BdG 
equations, the operators 7^ = "f_ E , and so these pairs of ±E 
energy states are not independent 2 . This is a consequence of 
the doubling of the number of degrees of freedom in ([Tp 6 . 

The BCS ground state is a state that is annihilated by 
1e\BCS) for E > 0, and by j E \BCS) for E < 0, which 
really are the same set of operators. Naively, one would be 
tempted to think of the ground state as a fermi sea of quasi- 
particle levels, all occupied below the Fermi level E < 0, 
and all unoccupied for levels above E > 0. However, as we 
see here, the "hole" r f_ E \BCS) and particle -y E \BCS) exci- 
tations (here E > 0) are in fact the same, and so the the BCS 
ground state should not be thought of as a filled Fermi sea 
of Bogoliubov quasiparticles. In the physical interpretation 
of the spectrum of the BdG equations, only the eigenvectors 
with energies E > should be understood as wavefunctions 
of physical excitations. 

Despite the subtleties of the physical interpretation of the 
BdG equations in this form, the purely mathematical analysis 
of them as an eigenvalue problem, is extremely useful in iden- 
tifying topologically protected zero modes. The BdG Hamil- 
tonian (the matrix-operator of |2|) has a spectrum with ±E 
energy pairs. For zero modes (E = 0), if some perturbation to 
the BdG matrix-operator were to cause the zero mode eigen- 
state to acquire nonzero energy, there would have to be an- 
other zero mode eigenstate acquiring the opposite nonzero en- 
ergy. The Atiyah-Singer index theorerrP3 tells us that the num- 
ber of solutions to an eigenvalue problem cannot be changed 
when deforming the operator continuously. Plainly put, new 
eigenstates cannot appear out of thin air - they must be defor- 
mations of eigenstates of the unperturbed system. Therefore, 
if the zero modes appear in pairs, they can in principle split. 
However, if the number of zero modes is odd, then at least one 
zero mode cannot be split by any physical perturbation (those 
that preserve the general form of the BCS Hamiltonian ([T|), 
and it is thus topologically robust. 

The most celebrated example of this, is the p x + ip y state, 
in which there is one zero mode bound to a (unit vorticity) 
vortex^ core, and so it is topologically protected. If the BdG 
matrix-operator spectrum has some degeneracy mandated by 
symmetries of the system, that are conserved by physical per- 
turbations, then the same degeneracy must hold. It then suf- 
fices for the zero modes to have an odd number modulo the 
minimal degeneracy mandated by the symmetries of the sys- 
tem, in order to be topologically protected - deviation from 
zero energy would split the zero modes in half, and the de- 
generacy protected by the symmetry would be violated. It is 




Phase winding 



FIG. 1: The punctured infinite plane geometry, with a magnetic flux 
threaded through the puncture (a vortex). 



however important to realize that if the deformation breaks the 
symmetries dictating the degeneracy of the entire spectrum, 
the zero modes may split. The zero modes are therefore pro- 
tected by the combination of symmetry and topological pro- 
tection. 



III. EQUIVALENCE OF VORTEX BOUND STATES AND 
EDGE STATES 



In this section we will argue that quite generally the vor- 
tex core bound zero mode states should be identified with 
zero mode edge states, in condensate phases. The connec- 
tion between the presence of zero modes at vortex cores and 
at sample edges has been previously examined in the co ntext 
of p x + ipy condensates of quadratic-dispersion fermionsOSD 

Consider the infinite plane with the order parameter ampli- 
tude non-zero and radially uniform only in the range r > L. 
Note that in the case of the p x + ip y condensate the order 
parameter has the forrrpJ {A, d x + id y }, and we refer to A 
as the condensate amplitude. The order parameter will in- 
clude a phase winding A = |A|e lm ^ (here and throughout 
this section, m is an integer), so that the region r < L at the 
disc center models a vortex core. A state bound to the vor- 
tex core will have a radial profile for the wavefunction that 
is decaying exponentially ~ e~ Ar in the superconducting re- 
gion r > L. In order for the wavefunction to be normalizable 
in the 7-^0 limit as well, the wavefunction amplitude must 
have the form \tp\ ~ r v with v > —1, in this limit. For 
v = —1, the wavefunction norm \tp\ 2 rdr ~ log f | e — >o di- 
verges logarithmically. We have therefore uncovered another 
boundary condition on the wavefunction - tpr must vanish at 
r = 0. Since bound state zero modes should not be sensitive 
to the detailed boundary condition inside the vortex coreM the 
r = point can then be mapped onto a hard wall of a small 
radius r = 5, on which tpr = 0, and therefore tp = on this 
wall. An even simpler model of the vortex is obtained when 
we identify 5 = L. The geometry then consists of a punctured 
infinite disc (see Fig. [TJ, with the condensate order parameter 
amplitude having uniform magnitude, and including a phase 
winding. 

The geometry of the punctured infinite plane can be contin- 
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FIG. 2: The semi-infinite cylinder geometry, topologically equivalent 
to the punctured infinite plane. There is a phase winding around the 
cylinder, corresponding to the phase winding of the vortex. 




FIG. 3: The semi-infinite plane geometry, topologically equivalent 
to the punctured infinite plane. 



states exist, in the presence of a vortex. 

It is instructive to examine how the momentum quantiza- 
tion evolves when mapping between the vortex and the edge 
geometries. Consider the semiclassical limit introduced in 
Ref.|24] where the quantization of angular/linear momenta can 
be inferred from a Bohr-Sommerfeld quantization rule of the 
classical orbits 



efoc 



A(n) ■ dn 



(3) 



Where here and in what follows I is integer, and the second 
integral on the RHS is the Berry phase traced by the (classi- 
cal) Nambu vector fi. Note that as opposed to Ref.|24j we do 
not transform to the London gauge. This is the result for a 
single band of quadratically-dispersing fermions, and we as- 
sume that a similar quantization rule will appear when taking 
a similar semi-classical limit of the more general problem. In 
particular 



2tt£H = <Z> p • <ix + topologically invariant terms . 



(4) 



Then in the process of deforming from the vortex to the edge, 
only the first term on the RHS changes. For a rotationally 
symmetric system, if the the angular momentum of the vortex 
core bound states is quantized Rp — h(£ + 7) (0 < 7 < 1), 
then the linear momentum for the edge states will be quantized 
as q = (J. + 7) ( where L = 2nR is as before, the system 
size in the direction parallel to the edge). 

Typically, the edge state energies E ~ q at low mo- 
menta/energy, and so we will need 7 = or integer angular 
momentum states in the vortex core, in order to support zero 
modes. 



uously deformed into a semi-infinite cylinder geometry (see 
Fig|2]l, with radius L. Taking L — > 00 then turns the geometry 
into the semi-infinite plane, Fig(3] The boundary condition at 
r = L corresponds now to a hard wall sample edge ^| wa n = 0. 
The radial exponentially decaying solution of the punctured 
infinite plane geometry will deform into an exponentially de- 
caying solution in the direction perpendicular to the edge (see 
Pig [3j, and the decay length will remain the same in both ge- 
ometries (in terms of the physical length scales in the system). 

The phase winding can be ignored in the formal limit 
L — * 00, but must be taken into account when considering 
a finite size of the plane. The smooth deformation we em- 
ployed to map between the vortex and an edge, will also map 
A = |A|e lm * -> A = \A\e im27r T. The phase winding in the 
order parameter will change the bound state wavefunctions 
qualitatively. When going from y to y + L we will pickup the 
requisite phase of 2nm. 

With the insight from the last section, it is now clear that 
modulo symmetries that are not related to position space the 
mapping we have described here, identifies the vortex core 
topological zero modes and topological zero mode edge states. 
Furthermore, we may conclude that showing the existence of 
one implies the existence of the other. Indeed, for the p x + ip y 
superconducting phase it is knowrPEED that zero mode edge 



IV. CONDENSATE PHASES 

In this section we briefly present the variety of condensate 
phases we will be examining in this manuscript. 

We begin by considering a simple model of fermions on 
the honeycomb lattice, either spinless or including spin. We 
include nearest neighbor hopping (of strength t) and fermion 
density-density interactions 

+ V f f f v al3 f f f 

where fj a are the bare fermionic operators. The indices i, j 
run over the sites of the honeycomb lattice, and the greek let- 
ters a, (3 =| , J. denote the spin indices. We point out that we 
neglect the gauge field in all our calculations. For spinless 
fermions, the indices a, (3 should be dropped. The interaction 
matrix is symmetric V^j 3 = V^ a , and is chosen such that 
V^ a = 0, so that /i indeed will be the Fermi energy. 

Throughout this manuscript we will assume we are in the 
weak interaction limit f > V«, so that BCS mean field theory 
is applicable. 
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FIG. 4: Conventions for the honeycomb lattice. The (blue) vectors 
marked di,2,3 represent the nearest neighbor vectors from sublattice 
1 to sublattice 2, and the (red) vectors a%,2 mark the Bravais lattice 
vectors. In our conventions, the nearest neighbor distance is set to 1. 



A. Order parameters 

An order parameter for off-diagonal long range order can be 
chosen as A"- 3 = — 2V i " >3 ((fiufjp)^) and must be antisym- 
metric A"? = —Aj". For spinless fermions, or spin triplet 
(A symmetric in spin indices) condensates, the order parame- 
ter is anti-symmetric in the lattice sites, and must break parity 
(the order parameter becomes anti-symmetric in swapping the 
i, j indices). A mean field BCS Hamiltonian can then be ob- 
tained 

T-LbCS = -t^2 fLfja + A* ^2 fjafja 
(ij)ct ja 



(6) 



In momentum space the BCS mean-field Hamiltonian reads 



Hues - 



a 



where f^ a (q) are the bare fermionic operators in momentum 
space, and we have introduced the matrix 



r(q) 



7 (q) 

7(q)* o 



(8) 



The indices /i, v = 1,2 denote the two triangular sublattices 
of the honeycomb lattice. Note that while we use /i for both 
a sublattice index, and for the chemical potential, it should be 
clear from context when /i is used for one or the other (specif- 
ically, whenever both appear in the same equation, the index 
/i is always a subscript). Furthermore, 7(q) = X)/=i e +lq df 
where di.2,3 are the three vectors from any given site in sub- 
lattice 1 to it's 3 nearest neighbors on sublattice 2 (see Fig. [4] 
for an illustration of our conventions). Finally, momentum is 
summed over the first Brillouin zone. 



For spinless fermions or spin-triplet condensates, the parity 
broken Ay implies A^(q) = — A"^(— q). We now turn to 
several condensate order parameters of interest. 

First, we mention the spin-singlet condensate phases intro- 
duced in Ref.|2l] which take the form A^ = Aijia^g. Here 
and elsewhere we will use the notation a x,y,z for the spin Pauli 
matrices (see table[lji. The function Ay is then symmetric, and 
Uchoa et al. 21 take Ay = £y Aq for an s-wave order parame- 
ter, and Ay = Ty |Ai, where Ty is the adjacency matrix for 
the honeycomb lattice (takes a value of 1 for nearest neighbor 
sites, and zero otherwise). The Ai order parameter mimics 
the structure of the tight binding kinetic energy term for the 
honeycomb lattice, Yluj) ■ ■ ■ = \ Tlij ^ij ■ ■ •> an( ^ m f act me 
Fourier transform of xTy is simply the matrix T(q) of ((8). 
As a result, near half filling just as the tight binding term has 
two Dirac nodes, so does the order parameter Ai. The Ai 
order parameter then has the approximate form of a p x + ip y 
order parameter. It may be a bit surprising to find a p x + ip y 
pairing function in a spin-singlet condensate, since it is then 
anti-symmetric under both spin exchange, and momentum in- 
version. However, the additional structure from the sublattice 
basis provides a third antisymmetric component of the pairing 
function, that keeps the overall anti-symmetry. It is important 
to emphasize at this point that this p x + ip y phase is gapped 
only when we are away from half filling (/1 7^ 0). 

Now we turn to spinless/spin-triplet order parameters. It 
is most convenient to write the spin-triplet order parameter 



in the form Aff = 



ia y a ■ di 



Here the (3-component) 



a/3 

vector is antisymmetric dy = — dji. Let us focus on a sin- 
gle component of the vector dy, to simplify our analysis, 
and also because this is equivalent to the spinless fermion 
case. Let us denote this single component as Ay, which is 
still anti-symmetric. The Fourier transform of this pairing 
function is a matrix A„„(q). The sublattice structure now 
mimics the behavior of the spin matrix structure, and can 
allow both sublattice spinor-singlet as well as triplet struc- 
tures. We introduce a new set of Pauli matrices if^^ in 
the 2-sublattice space (see table [E|, and now the spinless 
fermion pairing function can be written in complete gener- 

As men- 



( ? ) ality as A M „(q) 



A(q)+irf A (q) 

tioned earlier in this section, the order parameter for spinless 
fermions will satisfy A M!/ (q) = — A ufi (— q), implying that 
A (q) = A (-q) and A(q) = -A(-q). 

The simplest momentum structure in the order parameter 
A A11/ (q) would be just a function uniform in momentum space. 
A valid antisymmetric s-wave order parameter for spinless 

^A \ yA 
zA ) = ^ A °' 
However, attempting to transform this order parameter back 
to real space reveals that it break the honeycomb lattice sym- 
metries. Namely, the pairing function Ay is non-zero only 
for i, j on different sublattices, and in the same unit cell. The 
other nearest neighbors pairs of either i,j do not enjoy a pair- 
ing amplitude, and so discrete rotation symmetry is broken. 

For an order parameter that is linear in momentum (at least 
in a continuum limit), and maintains all the symmetries of 



fermion pairing is A A1!y (q) 
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the honeycomb lattice ( apart from the inversion symmetry, as 
mentioned above), we take A.y nonvanishing only on nearest- 
neighbor links. We choose for the directed links from one 
sublattice to the other the value Ajj = A2, and opposite for 
Aji = —A-2- In momentum space, the order parameter yields 



A„„ = A 2 



7 (q) 
-7(q)* 



(9) 



This pairing function can be conveniently rewritten as A M „ = 
A 2 77 z r(q). We note that this pairing function includes parts 
that are symmetric as well as anti-symmetric in the sublattice 
space. This pairing function is a directed version of the link- 
pairing order parameter in Ref. |2T| and similarly, when the 
fermions in the system are near half filling, ap x +ip y structure 
appears from the matrix T(q), and the order parameter has the 
approximate symmetry of a p x + ip y order parameter. 

With the spinless p x + ip y pairing, the bulk energy 
spectrum for the precise lattice model is found to be E 



± 



(i 2 + |A 2 | 2 ) | 7 | 2 + ri 2 ± t|7lVIA 2 | 2 (7 + cc) 2 + 4/i 



1/2 



and it has nodes at the points q = p ( 1, 



where 



and symmetry related 



2Vt 2 -|A 2 |= 

points. These points satify 7 real, and | 7 | 2 = pzn^p ■ 
Regardless of where the Fermi surface is, the nodal points are 

near the Fermi surface defined by | 7 | 2 = since |A| <C t. 
A superconducting phase with nodes very close to the Fermi 
surface allows for bulk states below the momentum space 
averaged gap. Our analysis in appendix|A]will indeed confirm 
this. 

In the following sections we will perform an exhaustive 
analysis of the 2 spin-singlet pairing phases to determine 
whether they allow topological zero modes. The analysis of 
the spinless phase we leave to appendix [Ajbecuase this phase 
has gapless bulk excitations. 



B. Continuum limit 

Analyzing the bound states in a vortex or an edge, is most 
easily done in a continuum limit of the lattice model. Here 
we will follow the conventions of Ref. 28 and work in the 
so-called "valley-isotropic" convention of the honeycomb lat- 
tice fermionic models near half filling. We demonstrate how 
this continuum limit is used on the BCS Hamiltonian 25 for the 
various condensate states we will consider here. 

The Dirac nodes are a pair of points where 7 (±Q) = (see 
Fig. [5] for illustration). In our conventions, the Dirac nodes 

are positioned at Q = f^gjO^ and — Q. The fermi oper- 
ators are expanded about these two nodeP, in the so-called 
"valley-isotropic" convention, and the two modes are identi- 
fied as right(R) and left(L) 



/ MQ (q+ Q) ~ nVv*ii(q) 

/>ta(q - Q) ~ ni^^aifq) 



(10) 




FIG. 5: First Brillouin zone of the honeycomb lattice, showing the 
two Dirac node positions ±Q (in red), and the main (x,y) axes indi- 
cated (in blue). 



Here II is some normalization, and we have used the y-Pauli 
matrix rf acting in the sublattice (/1, v) spinor space, that was 
introduced in the previous section. 

In addition to the Nambu, spin and sublattice spinor Pauli 
matrices, we now introduce a fourth set of Pauli matrices 
t x ' v ' z that act in the Dirac valley {A, B — R,L) spinor space 
(see table [I}. For clarity, from this point on, whenever it is 
convenient we will suppress indices which are being summed 
over. 

We expand the kinetic energy about the 2 Dirac nodes, us- 



ing 7 (q ± Q) 
obtain 



T2 (PxTWy) 



(p x ±iPyY 



to 



r(q±Q) w (rfq x ± r?q y )+\ [rf (q 2 x - q 2 ) T 2rfq x q y ] 

(11) 

Organizing the expressions using the various Pauli matrix sets 
we have introduced, and using the continuum Fourier trans- 
form 



«M(q) = J d 2 



(12) 



we derive a real space continuum version of the kinetic energy 



d 2 x^ 



\i — iv (ff ■ V) 



+- a t z (rf (d 2 x - d 2 ) - 2 v yd x d v ) 



(13) 



where ff = x rf + y rf 1 , and v = t|. 

At this point we observe that keeping only the linear deriva- 
tives, we obtain the celebrated Dirac operator, and at that level 
of approximation the Kinetic term has an SU (2) symmetry of 
valley spinor (r x ' y ' z ) rotations, in addition to the spin SU(2) 
symmetry. The last term is a quadratic correction to the Dirac 
operator, that explicitly breaks the Dirac spinor SU (2) invari- 
ance, reducing it to a £7(1) symmetry of rotations about t z . 
This correction, while often ignored, will be examined in the 
calculations we perform here. 

Next we turn to the pairing term in the BCS Hamiltonian. 
With the three pairing phases we outlined in the previous sec- 
tion, and using the various Pauli matrices we defined to com- 
pactify the expressions, the lattice Hamiltonian pairing term 
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C. Bogoliubov-de-Gennes equations 



Hr, 



*£{/(q)[^ (Ap + A^q)) 
q 

+A 2V z T( q )]f(- q ) + h.c.} 



(14) 

In the continuum limit, we approximate BCS off-diagonal 
terms in the following manner 

£/(p)M(p)/(-p) 
p 

« E [/(q + Q) M (<i + Q)/(-q - Q) 

q 

+ /(q-Q)M(q-Q)/(-q + Q)] ( 15 > 

~n 2 E [^(q)^(q + Q)<i*(-q) 

q 

+ ^(q)H)^M(q-Q)Vii(-q) 

where M(p) is the pairing function. In our case, the pairing 
function from ( fl4| ) is most conveniently cast as 



M(p) = ia y (A + A^) + A 2 r) z r 



(16) 



Using the lowest order in the expansion (jTTJ, the identity 
rf (ff* ■ q) rf = — (rf ■ q), and absorbing a factor of | into 
both Ai 2> we find 

M(q + Q)iif w [ier y (A - A^rf • q)) - A 2 rf ??• q] wj* 
-wfM(q-Q) « [i^(-Ao+A 1 (j7-q))-A 2 ^rf-q]^ 

(17) 

Combining these results, we find the pairing term in the 
Hamiltonian becomes in the continuum limit 



H 



~E{« 



+T x A 2 i] z (rf- q)7] y V(-q) + h.c.} ■ 

(18) 

In order to be able to take slowly spatially- varying order pa- 
rameter amplitudes (A 0j i i2 X we need to reorganize the pairing 
term in real space as 



H 



1 



pairing 



i<T v T y ri v Ao — ia v T y ~ {Ai, (77 • V)ry s } 



+ir^{A2,^(rf-V) ?7 n]V'+/i-c.} , 

(19) 

where all the operators if) are function of x, and {..,..} denotes 
the anti-commutator. 



With the final continuum forms of the kinetic ( pj) and pair- 
ing (JT9j terms of the BCS Hamiltonian, we can derive the 
BdG equations following the details of section [IT] The BdG 
equations for the phases we examine in this manuscript take 
the form 



wi: )-b r 



(20) 



where HedG = Ho +Hi + H2 +H3 +H4. Here the kinetic 
term is 



Hq 



[fi-iv (rf- V)] 






[fi + iv (rf -V)] 



= /mj z - iv (rfd x + rfuj z dy) 
= fiuj z — ivD , 



(21) 



where we have introduced D = {jfd x + i~i y u) z d y ). The sin- 
glet s-wave pairing term is 



Hi 



+iA r) y a y T y 
-iA*i] y a y T y 

rj u tr v T v \A \ur e e i * oU ' 

rfa y r y \A Q \ [w x cos(0 o ) +w tf sin(0o)] 



(22) 



where — 0o is the phase of — i Aq. The singlet p x + ip y pairing 
term is 



^cr y r^i{Ai,(77- V)t7»} 



--CT y rVr7 y {|A 1 |ea;p + ^ 1 ^,l)} 



(23) 



where —(j>i is the phase of — iAi. Rnally, the spinless p x +ip y 
pairing term is 



H 3 = 







iT-i{A|,^(rT-V)^} 



-TVw a! {|A 2 |e +i * 2I ' ; *,Z)} 



(24) 



where — 2 is the phase of iA 2 . Knally, the quadratic correc- 
tion to the kinetic term is 



H, = V (ufrf {dl - d 2 y ) 2 v y d x d y ) 



-t z D*Dw z ti x 
4 



(25) 



This concludes our derivation of the BdG continuum equa- 
tions, which we will now investigate in a variety of geome- 
tries, to determine whether topological zero modes appear. 
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D. Quantization rule 



SNS boundary 



In Sec. Ill we pointed out how the momentum quantiza- 
tion evolves when mapping between the vortex and edge ge- 
ometries. We assumed that the Bohr sommerfeld quantization 
rule takes on the form of Q. In this short subsection we will 
briefly deduce what the quantization rule is for the specific 
cases we consider in this section. 

In the effective continuum theory, the Dirac valley spinor 
degree of freedom is independent of the real space position 
degree of freedom. Therefore, the only significant difference 
between the BdG Hamiltonians we deal with here and those 
dealt with in Ref. |24l are the appearance of additional spinor 
structures - the spin, sublattice and Dirac valley spinors. Fol- 
lowing the semi-classical derivation of Ref. |24] we can use a 
coherent state representation of spin g not only for the Nambu 
spinor, but also for the 2 other spinors as well. We intro- 
duce classical unit vectors for each one of the spinors h for 
the Nambu spinor, h for the sublattice spinor, and t for the 
Dirac valley spinor. For the cases where our fermions have 
spin, there is also a spin degree of freedom, for which we use 
s for the coherent state representation of the spin. The path 
integral formulation will include a Berry phase for each one 
of the unit vectors, so 



j> p • dx + - j> Au(fi) ■ dh+ - j> A a {s) ■ ds 



+ \ <£ A v (h) ■ dh + ^ (f> Ar(i) ■ di . 



(26) 



for the spinless case, the spin Berry phase term will be absent. 

Finally, following Appendix C of Ref. [24] the Bohr- 
sommerfeld quantization rule will be of the form Sb — 
2ir (£ + j)H (the 7 part coming from the order parameter 
phase winding). We find therefore, that indeed in the cases we 
consider here, the Bohr sommerfeld quantization rule takes on 
the form of Q. 



V. SNS JUNCTIONS 

Perhaps the simplest indication of whether topological zero 
modes can exist in condensate systems is when considering 
SNS (superconducting-normal-superconducting) junctions^. 
In this section we will explore the spectrum of states bound 
to an SNS junction in the various condensate phases we men- 
tioned in the previous section. This will serve a first step to- 
ward determining in which of these phases zero modes may 
appear. 

The SNS junction is modeled as an infinite strip of width 
L in the continuum limit, where the pairing function vanishes 
(see Fig. |6j. On the two sides of the strip are condensate re- 
gions, with a uniform pairing function, with a relative U(l) 
phase (f>. For the SNS junction we have the pairing function 



A(x) 




(27) 




FIG. 6: SNS junction on the background of the honeycomb lattice. 
Our conventions are that the SNS boundaries are at an angle a with 
the armchair y-direction in our conventions. 



where x is the coordinate in the direction perpendicular to the 
SNS junction walls. 



A. Singlet s-wave condensate 

Combining the kinetic energy ( f2T| and s-wave pairing term 
(J22j» the BdG equation for the s-wave condensate takes the 
form 



T-LBdGi 1 =Eip 

= jJLUS* 



ivD 



(28) 



A (u x cos (0) + uj v sin (</>)) v y T v a y 



1>, 



where (j> is the U(l) phase of the order parameter, and to avoid 
clutter we have dropped the subscript from both the phase 
and order parameter magnitude. It is convenient to work in the 
London gauge - we use a unitary transformation O = e _l ^ w 
in order to rotate the phase 0^0. 

We are free to choose tp to be eigenstates of r v ,a v , such 
that r v a v — + tu. This 4-fold degeneracy applies to the entire 
quasiparticle energy spectrum. We note here that the explicit 
appearance of t v , a y is misleading, because the BCS Hamil- 
tonian in this phase is in fact SU(2) invariant for both the 
spin and Dirac spinors. It is perhaps more appropriate to write 
a v T v = —e T e a a product of the two totally antisymmetric 
2x2 tensors in the spin and Dirac spinor space. The SU(2) 
invariance in both these spinor spaces then becomes evident 
(these appear in the pairing terms of the Hamiltonian). Fur- 
ther assuming that the SNS junction is aligned with a armchair 
line of the honeycomb lattice (y-direction in our conventions - 
see Fig. |6j and that ip only varies in the x-direction, which we 
are allowed to assume in a y-infinite system, we find the BdG 
Hamiltonian reduces to 



Ti-BdG = 



fiuj z — ivrfdx 



Au} x rj y To~ 



(29) 



Apart from the expected (see section [TTJ symmetry 

uj x 7iBdG^ x = ~T~£*BdG it i s a l so eas y t0 snow that 
T] x oj z HBdG r l Xu;Z = T^-BdG- This relation is special for y- 
independent states (once y-variation is allowed this is no 
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longer a symmetry of the Hamiltonian). The additional sym- 
metry will yield a further double degeneracy of the spectrum, 
so we expect every energy level to be 8-fold degenerate. The 
eigenstates of d29} obey 



1 



cos(</>) = — 1, which means only when the two superconduct 
ing slabs have a ir phase difference does a zero mode eigen 
state appear. Furthermore, we can find the spectrum of low en 
ergy states - we use the limit E <C A to approximate Eq. ([36] 
as cos (^p) + cos(</>) ~ which then yields 



d x tp = AtP= —jf [fiuj z -E- Auj x rj y Ta] ip 



(30) 



The matrix A is then diagonalized, using some (x- 
independent) similarity transformation U, and we find it has 4 
eigenvalues 

U~ 1 -A-U = 

-Diagonal f- \J ' A 2 - E 2 + in, A 2 — E 2 ~ in, 
v V 

+ ^A 2 - E 2 + ifj,, +y/A 2 — E 2 — tfi 



(31) 



The diagonalizing transformation is 

a iE 



a + iE iE — a iE — a a + iE 

—a — iE iE — a a — iE a + iE 

A<tt — AuT Aerr — Act 

Act Act Aerr Act 



(32) 



where a = \/A 2 — E 2 . When considering energy levels well 
below the condensate gap E <C A we have VA 2 - E 2 > 0. 
Therefore, for x > we have the normalizable solution 

*{x) = C/.(e(-^ SOT+ ^)5a 1 ,e(-^ S ^^)5a 2 ,0,0) 7 

(33) 

For x < —L we have the normalizable solution 

^(x)=O-U-(0, 0, e(+^ ZF+i ") S ^ , e (+VACT-i M ) f 62 

(34) 

We must now solve for the wavefunction ip in the normal re- 
gion, and then match the wavefunction at the interfaces. For 
— L < x < the pairing function vanishes, in which case 
A = j^t] x (/xw z — E), and since all solutions are normaliz- 
able in this region, we can write 

^(x) = e Ax ^(x = 0) . (35) 

Now we match the wavefunction at x = 0, —L. This will yield 
a set of linear equations with the variables dip, b\,2- It can 
then be recast as a matrix equation B ■ (a\, 02, 61, &2) T = 0, 
and for a non-trivial solution, we require that Det(B) = 0. 
The equation for the determinant turns out to be 



4cos(^)£; 2 4VA 2 -E 2 sin( 2 Jf)E 



2 cos 



A 2 A 2 
j -2cos(0) = . 



(36) 



V v 



The condition Eq. ( poj l induces a quantization of the en- 
ergy values. In particular, we can now investigate whether 
zero modes are possible. With E = 0, Eq. d36| becomes 



v I 1 6 
Ek — 2tt [n + - - — 
2L \ 2 2t 



(37) 



We find the spectrum is evenly spaced, with a spacing 5E = 
2X2^- 

The zero mode wavefunctions can also be found (from the 
null space of the matrix B, with E = and <j> = ir). Including 
the t and a spinors we have been ignoring, \ T = (1, ir) T and 
Xa = (I, icr) T we find a total of 8 solutions 



1 



x < -L 
-L < x < , 
>. x > 

(38) 

where ^(0) = Af (err, —ijar, 1, 77) , 77 = ±1 and A/" is a 
normalization factor. Here, and elsewhere in this section, the 
column vector 7/^(0) has the entries (ui, 1*2, Ui, i>2) T , where 
1, 2 are the sublattice indices. These two solutions are inde- 
pendent (even at half filling \l = 0), and in fact are also eigen- 



r\ar 



It 



states of ri z ijj x tjj vc , T = <TTtp naT and r] x u} z ijj vaT 
is also evident that ui x ip* a . T — aTip^, h ^ a ^ T . 

As mentioned in section IV B the honeycomb lattice tight 
binding model in the continuum limit includes a quadratic 
derivative correction to the Dirac operator ( |2~5| ). Since the 
degeneracy in t = ±1 stems from the SU(2) valley spinor 
symmetry, which is explicitly broken (and reduced to U(l)) 
by ( p5j ), we should investigate whether this term splits the 8 
zero modes we have found. This will be our first step in ex- 
ploring whether the zero modes appear in the original lattice 
model. Since we are considering here y-independent states, 
the correction reduces to 



H, = -r 



(39) 



The correction commutes with rfu) z , and is spin SU(2) in- 
variant, so the quantum numbers r\ and a are conserved, so 
only r can mix and a quadruple degeneracy of every energy 
level will still hold. A simple calculation yields that all the 
matrix elements between the zero modes induced by the cor- 
rection vanish, and so this term does not split the zero modes, 
in this SNS geometry, with the armchair alignment. 

So far we have only considered a very particular alignment 
of the SNS junction walls - the y-direction in our conventions 
for the honeycomb lattice. Now we turn to investigate whether 
different orientations of the SNS junction behave different. 
Rotating the SNS junction counterclockwise by an angle a 
(see Fig. [6|, and assuming the eigenstates only vary in the 
direction perpendicular to the SNS junction walls, which we 
denote by x', the only term that changes in ( |28"] > in the London 
gauge is 

D(a) = [cos (a)rj x + sin (a W] d x , = V x e +laujZ ^ d x > . 

(40) 
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It is easy to show that the unitary transformation U(a) = 
e -tf w >) w jjj j-Qf^g a _^ mapping this problem directly 
onto the problem with the SNS junction parallel to the y- 
direction U(a)^HBdG(a)U(a) = TtBdG(& = 0). Further- 
more, the symmetry [H, u> z rf] = will simply be replaced by 
[H(a),U(a) cu z T] x U(ay\ =0. Thus, we conclude that the 
eigenvalue spectrum, at least when ignoring the quadratic cor- 
rection to the kinetic energy, is completely independent of the 
SNS junction orientation, and will remain 8-fold degenerate. 

The quadratic correction to the kinetic energy with the SNS 
junction walls rotated, takes the form 



H 4 = |r* (b{a)*b{a) 



u z r] x 



(41) 



Now we want to examine how the correction transforms under 
the unitary transformation that rotates a — ► in H.BdG(&)- 
We note first that all terms in H4 apart from the operator -rf 
remain invariant under this unitary transformation. The cor- 
rection becomes 



= -t z [cos (3a)r] x uj z - sin (3a)rj v } d 2 x , . 



The emergence of the 3a factors may seem a bit surprising, 
but in fact this is a consequence of the underlying 3-fold ro- 
tation symmetry of the honeycomb lattice - the splitting is the 
same if we rotate the SNS junction by 27r/3. The correction 
naturally reduces to Eq. ( |39l > when a = 0. In fact, the first 
term above is precisely the a = splitting matrix multiplied 
by the factor cos (3a). Using this fact, we need only compute 
the matrix elements of the second term when projected onto 
the subspace of zero modes. The result is 



{^ n 'a'T'\U HiU^ar) = - sin (3a)5o./ 0-77^ — 



A 2 



LA) 



(43) 



In general the eigenvalues of this matrix will be non-vanishing 
(except for the special pathological cases — ^ = n and a an 
integer multiple of 7r/3), thus splitting the zero mode energies. 

To conclude this subsection, we have shown that in the SNS 
junction geometry, the honeycomb Dirac dispersion s-wave 
condensate can support topological zero modes, when an odd 
phase winding is present ( the ir phase difference between the 
condensate slabs). However, these zero modes split when we 
take into account quadratic corrections to the kinetic energy, 
which are intrinsically present in the honeycomb lattice. Only 
in the case where the junction walls are aligned as armchair 
boundaries in the honeycomb lattice (a = 0), do the zero 
modes remain unsplit by the quadratic correction. 



B. Singlet p+ip condensate 

Now we turn to SNS junctions in the p x + ip y singlet phase. 
Our analysis will be very similar to that carried out in the pre- 
vious subsection for the s-wave phase, and as such we will de- 
scribe our calculations in much less detail. Here and through- 
out the remainder of our manuscript, we will assume /i > 
for all of the p x + ip y phases, since these are gapped only 
away from half filling, and since essentially the same result 
can be found for \i < (due to the honeycomb particle-hole 
symmetry). 

With a (piecewise) uniform pairing function, combining the 
kinetic energy pT) and the spin-singlet p x + ip y pairing term 
(|2"3"j) the BdG equation takes the form 



T-iBdai 1 = Eip 
Huj z -ivD-A (cj x cos (0) + u y sin (</>)) rf £>t V 



(44) 



where we have dropped the 1 subscript from both the order 
parameter phase and magnitude, to avoid clutter. As before 
we will work in the London gauge - the unitary transformation 
O = e"^^' will rotate the phase 0^0. 

Given what we have learned about the significance of the 
SNS junction orientation in the previous subsection, we begin 
by briefly addressing this point. As in the s-wave case, we 
assume the angle between the SNS boundaries and the y-axis 
is a, and consider eigenstates with spatial variation only in 
the direction perpendicular to the SNS boundaries. We use 
the same unitary transformation U (a) to rotate a — ► in the 
kinetic energy. The pairing function now includes D, which 
is also rotated to it's a = value, and all the other operators 
remain invariant. Choosing in addition eigenstates of r v , a y , 
the BdG Hamiltonian then reduces to 



T^-BdG = ^ Z - ivrj x d x + iAu; x r] z Tad :l 



(45) 



At this point we will note, that the BdG Hamiltonian we obtain 
here, just as its s-wave counterpart, has an extra symmetry 
[Hsdc rfuj z \ — 0. As a result, we expect the spectrum to be 
8-fold degenerate. 

Following the same procedure elaborated in the previous 
subsection, we find the energy quantization condition in the 
SNS junction to be 







M 2 )A 2 



(E 2 - fi 2 ) cos(0)A 2 

/2LE\ 



A 2 (£-/i) 2 
-((E 2 

+ 2ivE v / (E 2 - n 2 ) A 2 + v 2 E 2 sin 



2v 2 E 2 ) cos 



\ v J 
(2LE 



\ v 



For low energies E -C A this reduces to cos (</>) + cos 
and we obtain the low energy spectrum 



V7T 



2tt 



(46) 



' 2LE N 



(47) 
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identical to the spectrum we found for the s-wave phase, and 
including zero modes, only when </> = tt. 

Using the same conventions we used for the s-wave case, 
the zero mode wavefunctions we find, when taking <fi = tt, are 



e«(»+i>lA) 

1 

R ~»(l>-i7,A) 



x < -L 
-L < x < 

x > 



(48) 

where ^(0) = A/" (— icrr, jrycrr, 1, 77) and A/" is a normaliza- 
tion factor. Quite similarly to the s-wave zero modes, these 
solutions are eigenstates of r] x ui z ip naT — —ijip^ar and of 
rf u> v 'i/jr/ar — cTip^crr. From the form of the solutions it is 
also clear that u) (ipr/ar) — wt^>—ti,-u,-t- 

Finally we discuss the influence of the quadratic correction 
to the kinetic energy. Here the value of a is significant, and 
so we go straight to \A2) , and calculate the matrix elements in 
the zero mode subspace. We find that all the matrix elements 
vanish, and so zero modes are not split to first order in the 
correction. 



C. SNS junctions summary 

To conclude this section exploring the SNS junction geome- 
try, we recap the results of our calculations. For simplicity we 
have limited our discussion to wavefunctions uniform in the 
direction parallel to the walls. The spin-singlet s-wave phase 
supports zero modes only when a = (armchair boundary), 
and otherwise does not possess zero modes, with splitting due 
to the quadratic correction to the kinetic energy. The spin- 
singlet p x + ip y phase supports zero modes (to first order in 
perturbation theory in the quadratic correction). 



VI. EDGE STATES 

In this section we will investigate the edge state spectrum of 
the various phases we are exploring in this article. For conve- 
nience we will consider an edge where the honeycomb lattice 
abruptly ends, and assume the pairing function is uniform in 
space. We expect the bound states with low energy to appear 
with low momentum in the direction parallel to the edge (in 
the lattice model), and because of this one needs to be some- 
what cautious when thinking about the continuum limit. For 
the armchair edge of the honeycomb lattice, the Dirac point 
momenta are perpendicular to the boundary, and so low trans- 
verse momentum can be well described even in the continuum 
limit. For a zigzag edge, the Dirac momenta Q are paral- 
lel to the edge, and so small momentum in the lattice model 
p = Q + q~0 parallel to this edge, corresponds to momen- 
tum of order the Dirac momentum in the continuum theory 
q m — Q. Under these extreme conditions, the validity of the 
continuum limit approximation for the lattice model breaks 
down - the real momentum is quite far away from the Dirac 
point. We will therefore explore only the armchair edge in our 
present work (corresponding to a — in the previous section 
and shown in Fig.|6]l. 



A. Boundary conditions in the continuum limit of the 
honeycomb lattice 

Since we are taking a continuum limit of lattice models 
on the honeycomb, we must study with some care how the 
boundary conditions must be taken in the continuum limit. 

With our choice of the armchair edge, the boundary condi- 
tion of the lattice wavefunction is that it must vanish on some 
line. The eigenstates of the system are in general the Bogoli- 
ubov quasiparticles, with creation operators 



7' 



^[fi, (r)/t tt (r)+v(r)/^(r)] , (49) 



and the boundary condition corresponds to u = v = at 
the system edge. We note here that our description applies to 
non-condensate system as well, by simply taking v = ev- 
erywhere, in which case the Bogoliubov quasiparticles simply 
become modes of the fermi gas. 

In the continuum limit we employ here, the Bogoliubov 
quasiparticles are 



fj,aA 



(50) 

Changing to momentum space in both cases, and using the 
transformation Eq. ( fTO) , we identify 



^a!fl(q) = iVa(q + Q)n 

u mL {q) = -iiil v u va ((\ - Q)n 
?Vai?.(q) = *v«(q + Q)n 

V*L(q) = +w?jL««*(q- Q)H . 



(51) 



Using these relations we can identify the continuum limit ap- 
proximation of u, v. We find 



u m {r) « II [u^ aR e +lQ - r - 



i-nl„u vaL e 4Qr ] 



(52) 



Next we will use this continuum approximation of the lattice 
wavefunction to explore how the boundary conditions trans- 
late in the continuum limit. 

In our convention, the armchair edge can occur at the line 
x = 0, and since we have Q = Qx, the lattice wavefunction 
boundary condition translates into the condition 



Ur 
Vr ■ 



—irfuL 
+iifv L . 



(53) 



for the continuum wavefunction, at x = 0. Using the various 
Pauli matrix sets we have defined earlier in this manuscript, 
we can reorganize these conditions into the simple form 



U! z rj y T y il> — ip 



(54) 



where ip — (u, v) ( the continuum limit wavefunction). In 
what follows, we will assume the lattice occupies the x > 
semi-infinite plain, and use the boundary condition we have 



12 



derived here. It is worth while noting that when dealing 
with non-condensate wavefunctions, where uj z tp = ±tp, the 
boundary condition we have derived here simply reduces to 
the previously derived armchair boundary condition® in the 
conventions of Ref. W\\rj v t v tb = ±ip. 



B. Singlet s-wave condensate 

We now turn to explore the edge states in the honeycomb 
spin-singlet condensate s-wave phase. Starting from the BdG 
equations Eq. ( |28] i, we choose the order parameter phase (p = 
0, and solutions that are eigenstates of <t v i/j = atp and r y ip = 
rip, in which case 

[hoj z - iv (r] x d x + u z rfd v ) - uj x rf>ToA\ tp = Etp . (55) 

In the semi-infinite geometry, the system is still translationally 
invariant in the y-direction, so we choose solutions of the form 
tp(x, y) = e lqy tp(x). The BdG equations then become 

[p,uj z - ivrfd x + vquj z rf - uj x t] v tctA} tp(x) = Etp(x) . 

(56) 

The boundary condition (|54j> then requires w z r\ v tp(x — 0) = 

Ttp(x = 0). 

The set of coupled ODEs can be solved in a manner very 
similar to the way we solved for the x > region of the SNS 
junction, we cast the equations in the form d x ip — Atp, with 
the matrix A being independent of x. We diagonalize the ma- 
trix A with a similarity transformation that is x-independent, 
and then keep those solutions that are exponentially decaying 
in x > 0. In contrast to the SNS junction case, here we allow 
for a transverse momentum, and for this reason the calcula- 
tions are somewhat more involved. These solutions can be 
written as 



/ ri(A 2 +ir,(B+iEri)(E-^))aT \ 
2A(B-ir]fj.) 
(B+iEri)(TT(qv+F v ) 
2A{B-iri[i) 
ri(qv+F n ) 
2(B-i w ) 

U J 



where we have introduced B 

\J (B — irjfi) 2 + (qv) 2 , and 77 



= VA 2 -E 2 and 



(57) 
F„ = 



ir]fi) 2 + (qv) 2 , and 77 = ±1. The components of 
the 4-vector above correspond to the wavefunction amplitudes 
(til, u 2, V2) T , where 1,2 are the two sublattice indices. 
Throughout this section all 4-component vectors will follow 
this convention. 

Note that for small energy B < Awe have B « A, and 
then further assuming that the momentum is small qv <C A 
yields F n « (A — 277/1), which gives the decay length scale 
we found for the SNS junctions (as well as for the vortex core 
case, as we will see in the next section). At this level, before 
we impose the edge boundary conditions, we find we have 8 
solutions per energy E and transverse momentum q. 

The general spin-eigenvalue solution that decays exponen- 
tially in a; > is tp = tp+ia+i + ip-id-i. Now we will 
impose the armchair wall boundary conditions (|54j to this so- 
lution. The boundary conditions gives 2 linearly independent 



equations in the variables a n . These equations can be cast in 
matrix form, and for a non-trivial solution a n ^ to exist, the 
determinant of the matrix must vanish. The resultant equation 
for the determinant, is a quantization condition for the ener- 
gies E. The precise form of this quantization rule is 



=irB 3 + E11 {F-x - F x ) 
-iB(q 2 rv 2 + q(r(F- 1 + F x ) -2E)v 

-E (F_i + Fx) + t (F-xFi - fi 2 ) 



(58) 



Considering low energies E <C A, and small momentum 
qv <C A, we can approximate the quantization condition to 
{(E - qvr)A 2 + Efi 2 ) w which yields 



E 



qvA 2 



(59) 



We find that zero modes exist for q = 0. Next we will ob- 
tain the zero mode wavefunctions, by taking q = 0,E = 0, 
we recover the amplitudes (a + i, a^i) = (i — r, r + i). The 
complete wavefunctions of the zero modes are 



(60) 



where ipo = (<jt, —ar 1 1,1) . We find a total of 4 zero modes 
(a,r = ±l). 

The particle-hole relation of the BCS Hamiltonians that 
given an eigenstate ipE with energy E, uj x i\)* e is also an eigen- 
state with energy —E, when applied to the 4 zero modes we 
find here, gives 4 states that are orthogonal to the zero modes 
we found. This surprising result is understood when con- 
sidering how the boundary condition behaves. Starting from 
oj z rj v T y tp = ip, we want to know what boundary condition is 
satisfied by -0 = w 1 ^*- It can be easily shown that the bound- 
ary condition is u> z rj v T y ip = —ip. This result shows us that 
the tp states are precisely the ones discarded by the boundary 
condition in this case! Therefore, the only states satisfying 
the boundary conditions are the 4 zero modes we have found 
above. Finally, it is amusing to mention another consequence 
of these boundary conditions - the superposition yielding a 
Majorana fermion ip + tp, cannot be taken here! Therefore, 
despite the existence of zero modes, they cannot form Majo- 
rana fermion states. 

The edge states energy spectrum at low momentum q is lin- 
ear in the momentum, and we now proceed to briefly calculate 
the approximate edge state wavefunctions for these low ener- 
gies. The boundary conditions, cast as linear equations in the 
coefficients a n , can be linearized in energy and momentum. 
In this case the energy quantization condition we derive is 



E 



Aqrv (qv + 2A) 



A (qv + 2A) + 2/j 2 - qv^T 



(61) 



with co = ±1. Linearizing in momentum q, this result reduces 
to ( |59) , The solution for a v we find with this value of E is 
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a+i = a*_ x = (A - ifi)(qv + (A + ifi)(%T + 1)). With these 
coefficients, linearizing everything in momentum q yields 



1pa,r = X<t®Xt 



-+iqy 



• x ^a. 



where 



/ (A 2 + i(A + iE)(E - fi)) or 
-(A + iE)(qv + A - ifj,)ar 
A(qv + A — i[i) 

V A(A - ip) 



(62) 



(63) 



We find that every energy level has a 4-fold degeneracy, in- 
cluding the zero modes. This is the minimal expected degen- 
eracy, required by the SU (2) symmetries of both the spin and 
the valley spinor. Since the zero modes posses only this mini- 
mal degeneracy as well, it is topologically protected. 

Considering the quadratic correction to the kinetic energy 
term, since the zero modes only vary in the x-direction, the 
splitting term reduces to ([39]) A straightforward calculation 
then shows that all the matrix elements between the zero mode 
wavefunctions i/;® T vanish, and we find that there is no split- 
ting to first order. This is in agreement with subsection |V A| 
because the armchair edge corresponds to the angle a = in 
(|43|>, yielding no splitting. 



C. Singlet p+ip condensate 

In this subsection, we turn to explore the edge states in the 
honeycomb spin-singlet condensate p x + ip y phase. Starting 
from the BdG equations Eq. ( p4| ), we proceed with a calcula- 
tion that is only slightly different than that performed for the 
s-wave phase in the previous subsection. 

The system is translationally invariant in the y-direction, so 
we choose solutions of the form ip(x,y) = e >qv tp(x). Fur- 
thermore, we choose the order parameter phase <f> = 0, and 
solutions that are eigenstates of <r v ijj = atp and t v i/j = rip. 
All this yields 



pLLO 



ivD - Auj x r) y Dra ip = Eip , 



(64) 



where D — (r] x d x + oj z i] y iq). And the boundary condition 
( |54| i then requires u> z r] v tp(x = 0) = rip(x = 0). 

We find the solutions to this set of coupled ODEs in the 
same manner as in the previous subsection. We find those so- 
lutions that are exponentially decaying in x > 0, and linearize 
them in momentum and energy, expecting a low energy rela- 
tion E ~ q. We find the general (linearized) solution 



( i(ivE + Ar](E + i±))<7T 
(vErj + qA(v - iAri)ri - iA(E + [i))<7T \ (6:,s 
Arj(q(A + ivrj) + (i) 
V A M 

where r) = ±1. We take the general solution ip = ^ i> v a v 
and find which coefficients a„ will satisfy the armchair bound- 
ary conditions. The set of equations for a v can be cast in a 



matrix form, and for a non-trivial solution to exist, the matrix 
determinant must vanish. This yields the approximate quanti- 
zation condition for the low energy spectrum 



2Acr (A (v 2 + A- 

= 0, 

with the solutions 
A 



<f 



(Ev 2 + A 2 (E + 2/i)) q - 2vE^t) 



(66) 



E = — 



q 2 v 2 



2qA[i) _ A 2 rq 



qv A 



2[iTV + qA 2 



+ 0(q 2 ) 



(67) 

As expected, we indeed find a branch of low-energy states 
with energy linear in the transverse momentum, and zero 
modes for q = 0. Next, we find the approximate edge state 
wavefunctions. The solutions for the coefficients are 



q(A — r/iv) + fj, + ifir/T 



(68) 



The full edge state wavefunctions we find then, after lineariz- 
ing them with respect to the momentum q, are 



COS 



- sin 



Xa- ® Xr 



VXfJ, 



A 2 



VXfl 



A 2 t) 
-A 2 ) 



A 2 



- ,», +iqy 

v 2 +A 2 

1 a {q (tv 2 — Av - 
ia (<? {y 2 — Atv 
iv(qA + [i)t 

\v(qA + n) 

—a {qA 2 + Vfxr) 
ia(vfi + qA(2v + At)) 
—iv(—2qA — fi + qvr) 
v(qv — /xr) 



Vfl) 

vfir) 



(69) 



We find every low-energy state is 4-fold degenerate - for E ^ 
0, there is spin degeneracy, and a 2-fold degeneracy of the 
product rq. For the zero modes (q = 0) there is still a 4-fold 
degeneracy, in both spin and valley spinor degeneracy. As for 
the s-wave case, this is the minimal degeneracy of each energy 
level, and as such, the zero modes are topologically protected. 

Finally, we turn to examine what influence the quadratic 
correction to the kinetic energy has over the zero modes, since 
this perturbation breaks the valley spinor SU(2) symmetry. 
As in the s-wave case, the zero modes have no y-dependence, 
so the correction reduces to ([39]). A straightforward calcula- 
tion of the matrix elements between the different zero modes, 
yields these all vanish, and so there is no splitting from this 
correction, to first order. 



VII. VORTEX CORE ZERO MODE BOUND STATES 

In this section we will investigate whether zero mode bound 
states at vortex cores exist, in various phases. 



A. Singlet s-wave condensate 

The simple s-wave singlet-pairing condensate phase on the 
honeycomb lattice, has an eigenvalue spectrum determined by 
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the BdG equation 



The full wavefunction solutions we find are, for 77 = +1 



(70) 



where we refer to pi) and |22) for the full details of the ki- 
netic and pairing term. In this section we will find the zero- 
mode solutions in ( f70] > explicitly. As opposed to the calcu- 
lation of Ref. [T9J we allow for a non-zero chemical potential 
(corresponding to slight deviations from half filling), and we 
find exact solutions for the zero modes. 

As a first step, we choose solutions that are a v r y eigen- 
states, precisely as in the SNS and edge geometries, and 
as a result the energy spectrum will be at least 4-fold de- 
generate. We will model the vortex by assuming the form 
-iA = A(r)e^, with A(r) real. 

We begin by considering the half filling case (77 — 0). The 
BdG equation then become 



-w(rf-V) A(r)e^r? y crr 
A(r)e-' i V^ -w(?T-V) 



ip = 



(71) 



In this special case, one can obtain solutions that exist on only 
one of the two sublattices. this becomes obvious when mul- 
tiplying the equation set Hip = by r\ v on the left, resulting 
in 



-iv (—iri z d x + d y ) A(r)e l ^(7T 
A(r)e~ l ^(TT — iv (—ir] z d x — d y ) 



ip = 



(72) 

We can then choose to consider a solution on one of the sublat- 
tices 1, 2 in which case rf — > 77 = ±1. Writing the equations 
in polar coordinates 

-v V e+^ (d r + ^) A(r)e+*or \ 

(73) 

We observe that a for 77 = +1, we can try a solution ip = 
(u, v) where u, v are independent of the angle <p, and for 77 = 
— 1, we can try a solution of the form ip = ^(e l ^u, e - *^u). 
For both choices, the equations reduce to 



— vd r A(r)riaT 
A(r)r]<jT —vd r 



= 



(74) 



Rewriting these equations using the Nambu pauli matrices 

[-vd r + A(r)r]<7TLo x ] {u, v) T = , (75) 
it becomes clear the solutions are uj x eigenstates, 



u 
v 



= 1^1 e T,,TTU> " fo dr ' A ( r ') 



(76) 



Already at this point, before taking into account normaliz- 
ability, we see that as many as 16 zero mode solution exist, 
parametrized by 77, a, t, oj = ±1. Assuming that A(r) is pos- 
itive at r — > 00, only half of these 16 solutions are expo- 
nentially decaying rjaroj = — 1, and thus normalizable in an 
infinite system. The number of zero modes is then reduced 
to 8. In a finite system the exponentially growing solutions 
correspond to edge states. 



' *j ^ \r \, I ° I e +^^/o r ^M , (77) 



and for 77 = — 1 



^2 = Xr ® Xo 



( \ \ 





e -»™i fo dr ' A ( r ') . (78) 



Here, as in the previous section, the 4-component vectors cor- 
respond to (ui, u 2 , vi, V2) T (1,2 are the sublattice indices). 
We will follow this convention in the remainder of this sec- 
tion. With the condition of normalizability, the exponentials 
must take on the decaying form e - " Ja dr A ^ r \ and we must 
have ui = —crrr). The solutions then become 



(79) 



where ^? = (1,0, err, 0) T and ip% = £ (0, e^, 0, crre"^) 7 '. 

Next we turn to normalizability in the r — > limit. This is 
determined by whether the integral \ ip \ 2 rdr diverges. The 
solution %p\, is clearly normalizable in this region, while ip2 
clearly is not. This leaves us with 4 zero modes, rather than 
8 as we found for the SNS junction, where no analog of the 
r — > normalizability condition appears. 

Now we turn to the case away from half filling. The BdG 
equations in this case read 



H-iv(rf-V) A(r)e l V<>-T 
A(r)e^crT -fi- iv(rf* • V) 



tp = 0. 



(80) 



We can eliminate the <p dependence from the problem by 
choosing the exact same form as in the half filling case 



V<r,0) = ( Ul (r),e' i %(r),7; 1 (r),e-^7; 2 (r)) J 



(81) 



The reduced ODEs then involve only the radial coordinate. 
At this point it is useful, to make the educated guess 



/ u x (r) \ 



ip(r, 4>) 



oi<t> 



u 2 {r) 
vi(r) 
\ e-**«2(r) J 



,-i/J-dr'A(r') 



(82) 



inspired by the form of the solution for the SNS junction. 
Plugging this form into the ODEs, does not remove the or- 
der parameter from them completely. However, choosing 
vi(r) = —arui(r) and V2{r) — <jtui(j) in addition, does 
remove the order parameter. The reduced ODEs, involving 
only tii 2 then rea d 



777*2(7") — ivu x (r) = 
77; 77,2(7") 



77771 (r) 



ivu' 2 {r) = 



(83) 
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Extracting U2 from first equation, and plugging it into the sec- 
ond yields a single ODE for u\ 



u[(r)v 2 u'l{r)v 2 



rfi 



/' 



+ (j,iti(r) = 



(84) 



The solutions are Bessel functions J ( 7 ^),Y ( 7 ^), and U2 is 

ivu 1 (r) 

we obtain are then 



obtained from U2(r) 



The two zero mode solutions 



= Xt 



Xt 



■i/J-dr'A(r') 



-°tJ (v?) 



(lb. 

A v. 



l -)J 



(85) 

and the second solution simply has all the Bessel functions of 
the 1st kind replaced with Bessel functions of the 2nd kind, 
with the same parameters. Only the Bessel functions of the 
1st kind is normalizable in r — ► 0, or alternatively (specif- 
ically Jq° Y\ ( ^ ) rdr diverges ), if we impose a boundary 
condition at some small r = a, we will pick out some combi- 
nation of the two Bessel function kinds. With the a, r degen- 
eracy we end up with 4 zero modes. It is easy to verify that the 
BCS particle-hole relation yields lo x (V'i 2) = ~ <JT ' l Pi2 , ~ r 
producing no new zero modes beyond the 4 mandated by the 
system symmetries. 

It is noteworthy that the Bessel function, far from the vortex 
core ^ ^> 1, has an oscillatory nature, with a length scale -, 
precisely as in the zero modes we find for the SNS and edge 
states. 



The reduced ODEs then involve only the radial coordinate, 
and can be cast in the form d r ijj — Aip where 



.4 



1 

27 



1 



A 2 



AfiaTUj y r] 

A 

27 



IVflT] U! 
,,2 



V(JTLO X r) X 



,v 
2r 



(88) 



We first find the asymptotic (r — > 00) solutions to the ODE 
system. Neglecting all the - terms, we find 



A = - 



1 

■j 2 + A 



2 [A[M7TUJ y r] Z 



IVflT) UJ 



(89) 



We can diagonalize the asymptotic form of A with the unitary 
transformation 



= 



(90) 



yielding 

O^AO = 



A 2 



Diagonal(— iv — Act, +iv — Act, - 



-iv 



Aijt, +iv + Act) 
(91) 



Using this unitary transformation on the full matrix A, we find 
a block-diagonal form 



O^AO = 



A aT 
A- aT 



(92) 



where 



B. Singlet p+ip condensate 

We turn now to the simple p x + ip y singlet-pairing conden- 
sate phase. The vortex core eigenvalue spectrum is determined 
by the BdG equation 



(Ho + H 2 ) $ = EiP , 



(86) 



where ( f2"Tj ) and |23| contain the full details of the kinetic and 
pairing term. 

As in the s-wave case, we choose solutions that are a v T v 
eigenstates, precisely as in the SNS and edge geometries, and 
as a result the energy spectrum will be at least 4-fold de- 
generate. We will model the vortex by assuming the form 
iAi = A(r)e +l< ^, with A(r) real (different from our conven- 
tions in earlier sections so that we can use the same ansatz 
for the polar angle dependence as for the s-wave case). Also, 
since it will prove convenient, we will assume that the order 
parameter radial profile is piecewise uniform - vanishing in- 
side the vortex core, and constant outside it. 

With the insight gained in the previous subsection, we find 
the 4> dependence can be eliminated from the zero-mode prob- 
lem by choosing the wavefunction form 



i>(r,cf>) = (« 1 (r) J e i *u 2 (r),«i(r) > e- i *«2(r))' J 



(87) 



•It 



-A/iCTT 



-IVfl 

vL _ 

2r 



iAarv 
2r 



2r 



2r 



(93) 



Here the terms outside the matrix are implicitly multiplied by 
2 by 2 identity matrices. 

We now turn to solve the reduced ODE system A aT ^ = 
d r £, where £ = (fi(r), f 2 (r)) T . The precise solution will 
prove cumbersome to work with, and so we will start with 
an approximate solution that will reveal all the features of the 



solutions we need to find. First we write £ 



7? c 



We note at this point that for the solution to be normalizable 
at t — ► 00, we must have /icr > 0, if however this is not the 
case, then we simply choose the solution for the A- aT sector, 
in O^AO. The ODE for C is then 



1 



A 2 




2r 



2r 

ivfj, 



C- (94) 



Now we assume that v ^S> A, and consider the small r limit, 
so that we can approximate 



1 

27 



1 

1 



c 



(95) 
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for which the solutions are 

C, = (l,r,) T r-"/ 2 . (96) 

In the limit r — > 0, only the T) = +1 solution is normalizable. 
We therefore find approximate solutions, that have only the 
a, t 4-fold degeneracy. 

We now turn to briefly make connection with the precise 
solutions for the zero modes. From the equation d r f\ = . . ., 
we extract / 2 (r), and substitute it in the other equation. This 
yields a single 2nd order ODE 

(A 2 + 8r^(TTA + irfi(iv + r/i)) /i(r) 

+ 8r (v 2 + A 2 + rAfiar) f[ (r) (97) 

+ 4r 2 (v 2 + A 2 ) fi(r) = . 



rfj.(iv + Aa-T) - ( — _ 

Next, we write fi(r) = g(r)e ^T^ 2 r 2 VV» 2 +a 2 
and furthermore replace the radial variable with r 

iz(v 2 +A 2 ) 
2vJl 



The ODE for g(r) is then 

z + l) g'{z)+zg"(z) = 0, 



2Vv 2 + A 2 '' \Vv 2 + A 2 

(98) 

which is the confluent hypergeometric ODE. The solutions are 

g(z) = c x M (a, 1 + 2a, z) + c 2 z~ 2a M (-a, 1 - 2a, z) , (99) 
where ci 9 are free coefficients, a = - , " ., < I (and 

' 2\/v 2 +A 2 — ^ 

also a > 0), and M(a,b,z) = iFi(a,b,z) is the con- 
fluent hypergeometric function of the first kind (or Kum- 
mer function). With the complex valued variable z ~ ir, 
the solutions are well-behaved at r — > 00. At small r, to 
lowest order M(a, b,z) 1 + O(z), and the radial part 



of the wavefunction behaves like /i(r) 



„— 1/2 a 



,,-1/2 a 



Cl 



Cir 



a-l/2 + c r -o-l/2 



c 2 r~ a ~ 1//2 . When a < | this solution is normalizable. How- 
ever, in the r — > limit the order parameter must vanish, so 
we must take A — 0, in which case a = 1/2, and causes a 
logarithmic divergence when we try to normalize it. There- 
fore, including the r — * point, we must set c 2 = 0, and we 
are indeed left with only 4 solutions for zero modes. 

C. Quadratic correction for the vortex core zero modes 

In both cases of the s-wave and p x + ip y spin singlet phases, 
the eigenstates are angular momentum eigenstates as well, and 
have the general separable wavefunction form 

/ ui(r) \ 



VvO, 4>) = X* ® Xr ® e l 



vi(r) 



(100) 



We take the quadratic correction ( |25] l in polar coordinates and 
find that it takes this wavefunction into the form 



H4ipi(r, ^>)=Xa® X-r <8 e 



it<t> 



-2i4> 



Mr) 
' 4 9i{r) 
P 92(r) J 



(101) 



Then trying to take the product (<ipt> \H4\1IJe), it suffices to con- 
sider the angular dependency 



( e^hir) \ 
5 -w*/ 2 (r) 

V e 2 ^. 92 (r) / 



2tt 



vl gi e-** + v*g 2 e 3 ^ 



(102) 



from the polar phase integration we conclude that nonzero ma- 
trix elements exist only between states with angular momen- 
tum £ differing by ±3. In particular, to first order, there is no 
correction, and the zero modes persist to this order. 

To conclude this section, we point out that the pure Dirac 
theory approximating the honeycomb lattice allows for topo- 
logical zero modes to appear bound to vortices in both spinfull 
condensate phases. 



VIII. NUMERICS ON THE HONEYCOMB LATTICE 
MODEL 

In the previous sections we found zero modes exist in both 
the s-wave and p x + ip y spin-singlet states in the continuum 
Dirac approximation for the honeycomb lattice. We also tried 
to ascertain whether the zero modes exist beyond the approx- 
imate Dirac theory for the honeycomb lattice, by taking into 
account the quadratic correction to the kinetic energy Eq. ( |25| . 
We calculated whether this correction splits the zero modes at 
first order in perturbation theory. With the exception of one 
case, we always found that to first order, no splitting occurs. 
The exception is the s-wave phase in the SNS geometry with 



a^0 (Section V A 1, where we found the correction does give 



a splitting to first order in the quadratic correction. In contrast, 
in the edge state and vortex cases, no such splitting occurred 
at first order. While the SNS splitting does vanish when the 
junction boundaries are of the armchair edge type, consistent 
with the edge state result, the collection of these results is in- 
conclusive as to whether the zero modes really do appear in 
the lattice model, and not just in the idealized approximate 
Dirac theory. To answer this question definitively, we have 
performed numerical calculations (exact diagonalization) on 
the precise lattice models for all phases where we suspect zero 
modes occur. 

We consider the vortex state case of the two spin-singlet 
phases at precisely half filling. First, we constructed lattice 
patches of square, rectangular, circular and elliptic shapes (see 
Fig. 1 1 for illustrations), of various sizes. We then diagonal- 



ized the matrices describing the lattice model |6]l on these lat- 
tice patches, with the 2 spin-singlet order parameters includ- 
ing unit-vorticity vortices at their centers (one representative 
example is shown in Pig J7]>. In all cases, we find the lowest 
energy eigenvalues Eq, and compare them with the de Gennes 
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FIG. 7: (color online) Circular lattice patch in the vortex state. The 
(black) arrows represent the phase at a lattice point as the angle the 
arrow makes with the x-axis in the picture. The (red) dots are the 
honeycomb lattice sites, and (blue) lines are the nearest neighbor 
links. Here, for illustration purposes, the radius of the vortex core 
is taken to be 2 (the radius of the lattice patch is 6, where the nearest 
neighbor distance is 1/ \/3). Following the arrows it can be verified 
that the vortex indeed has a unit vorticity. 
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TABLE II: Lattice patch sizes 



energy scale -g^, the energy scale one expects for vortex core 
bound state 20 . For the s-wave state, the gap energy E g = A, 
while for the spin-singlet p x + ip y state, the gap scales with 
the chemical potential 21 E g ~ /i. The sizes of the various 
lattice patch geometries we take is detailed in table [II] where 
the nearest neighbor distance is 1/ v3 (the primitive Bravais 
lattice vectors are then of length 1). The scaling of the lowest 
energy with the finite system size is described in Fig.[8]for the 
s-wave state, and in Fig. 10 for the spin-singlet p x + ip y state. 



It is clear from all scaling plots that the lowest energy is of 
order of the de Gennes energy, and that this energy does not 
significantly decrease with growing system size. This would 
indicate that there exist no zero modes in these phases, despite 
the results from the continuum Dirac theory. 

In addition, we plot the spatial density of the lowest energy 
quasiparticle density on the lattice patch, in a number of rep- 
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FIG. 8: Scaling of the lowest energy with lattice patch size. We plot 
the ratio of the lowest energy to the de Gennes energy for all sizes 
we take. The square lattice patcl |(a)| is of size L x L, the circular 
lattice patc r|(b)| is of radius L, the rectangular lattice patc r|(c)| is of 
dimensions L X |L, and the elliptic lattice patct (d) is of main axis' 
L and 2L. The corresponding number of sites for each lattice patch is 
detailed in table [n] It is clear from the graphs that the lowest energy 
is always of the order of the de Gennes energy, and remains roughly 
unchanged when we increase the system size. 



resentative cases in Fig. 1 1 and Fig. 12 in order to verify that 
these indeed are vortex core bound states. Specifically, at each 
lattice site j we plot a dot with its color signifying the relative 



magnitude of 



(the values are normalized to run 



between and 1, and the color is varied linearly with with 
this value). From the plots it is clear that these lowest energy 
states are indeed vortex core bound states. 

In all the cases described in Fig. [8] and Fig. 10 the vortex 
core size was taken to be 0, forcing the vortex phase winding 
to occur over a distance that is comparable to the lattice length 
scale. This fact is what invalidates the Dirac continuum the- 
ory - the order parameter in these cases is not a slowly varying 
function on the lattice scale, near the vortex core. Following 
this last observation, we also calculated the energy spectrum 
for a series of different vortex core sizes, ranging from 0.2 to 
5.8 in increments of 0.4, while keeping fixed the overall sys- 
tem size (circular lattice patch of radius 12.0). The results are 
plotted in Fig. [9] and it is clear from them that the energy split- 
ting decreases rapidly with the vortex core size. The highest 
energy we find (at the smallest radius) is 0.88 times the de 
Gennes scale, and the smallest energy scale we find is 0.0016 
(at a radius of 5.0). However, the de Gennes energy scale in 
terms of the correlation length of a superconductor is actu- 
ally 20 . From this we expect that the lowest energy should 



change with the vortex core radius R as Eq ~ j^. We fit the 
raw data in Fig.|9]to a curve — 0.0007+0. 17/R, also displayed 
in Fig. [9] When the vortex core size is bigger, the phase wind- 
ing occurs over a larger distance, and the approximation of a 
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FIG. 9: Ratio of lowest energy scale to de Gennes energy scale, with 
varying vortex core size. The vortex core sizes range from 0.2 to 6.0 
in increments of 0.4, and the lattice patch is circular with a radius of 
6.0. The nearest neighbor distance is 1 /v3- The raw data is denoted 
by the continuous (red) curve, and the fit is denoted by a dashed 
(blue) curve. 



slowly varying order parameter improves, but the lowest en- 
ergy is still of the de Gennes scale. We find therefore that the 
zero modes are split, and correspond to the de Gennes bound 
state spectrum. 

For a vortex one would expect the core size to be of the 
order of the correlation length in the superconducting state. 
The parameters we use are A = 0.5, t = 1.0 the lattice con- 
stant a — 1.0 and working in units where h = 1, we have 
a correlation length that is of the order of the lattice constant 
£ ~ TK ~ lis ~ 1 - In Fig. [5] we can see that for this core 
size, the lowest energy is between 5 and 10 times smaller than 
the de Gennes scale. 

In conclusion, the numerics we have done show that the 
zero modes appearing in the continuum Dirac theory are split 
in the lattice model, but that the bound state energy can be 
significantly smaller than the de Gennes scale. 



IX. EXPERIMENTAL REALIZATIONS - BOSE FERMI 
MIXTURES 

Superconductivity, which appears in many conventional 
fermionic systems, does not seem to occur intrinsically in 
graphene, the most readily available realization of the honey- 
comb tight binding model, but can be induced via the proxim- 
ity effect close to another superconducting materiaP^El The 
peculiar band structure of the honeycomb lattice is not how- 
ever limited to graphene - it is just one material realization 
(other possibilities may include thin films of quasi-2D honey- 
comb layered superconductorsSEEO) Another possibility is a 
cold fermion gas trapped in an optical lattice, with a fermionic 
atom density corresponding to about half filling. Since the 
atoms are electrically neutral, the interactions between them 
are to a good approximation simple collisions, corresponding 
to an on-site interaction in a lattice model. Superconductivity 
requires some attraction between fermions, so in order to have 



l'-IE F 

1.4c 

1.2 

1.0! 
0.8 
0.6 ! 
0.4 

0.2 : 



15 



20 

(a) 



25 



L 







1.4 




1.2 




1.0 




0.8 




0.6 




0.4 




0.2 




L 0.0 





10 12 14 16 

(b) 




L 





l' 2 /E F 


1.4 




1.2 




1.0 




O.S 




0.6 




0.4 




0.2 




L 0.0 





10 12 



(C) 



(d) 



FIG. 10: Scaling of the lowest energy with lattice patch size. We plot 
the ratio of the lowest energy to the de Gennes energy for all sizes we 
take, for the spin-singlet p x + ip y state. In this case, the gap energy 
scales like the chemical potential fi, so we take as the de Gennes 
energy scale fj, 2 /Ep. The square lattice patct (a) is of size Lx L, the 
circular lattice patc r|(b)| is of radius L, the rectangular lattice patc r|(c)| 
is of dimensions L x |L, and the elliptic lattice patcl (d) is of main 
axis' L and 2L. The corresponding number of sites for each lattice 
patch is detailed in table|n] It is clear from the graphs that the lowest 
energy is always of the order of the de Gennes energy, and remains 
roughly unchanged when we increase the system size. 



any hope of realizing such a phase, one needs to cause attrac- 
tive interactions between the fermions. In solids, phonons pro- 
vide this mechanism by inducing an attractive interaction be- 
tween electrons. In cold atom gases, this role can be assumed 
by adding bosonic atoms; sound modes in a Bose-Einstein 
condensate mimic phonons in a solid. The virtues of the cold 
atom realization do not end in simply making the supercon- 
ducting state feasible, but also provide a great deal of control 
over many parameters. 

Motivated by the reasoning discussed in the previous para- 
graph, we will now consider a model of a Bose-Fermi mix- 
ture on the honeycomb lattice model. We consider only on- 
site interactions, since usually one has to work quite hard to 
make longer range interactions appreciable compared to them 
in cold atom systems. Our Hamiltonian therefore reads 

n =- t Yl /U«+AEi^ 

(ij},ot j,ct 

-wJ2 b i b i +3w 12 b fa 

(ij) 3 

' " ' ' Pa 



u » E - po] + u f f E (/Ut) (fjjji 

3 3 

^EE(/U«) 0>0 . 



(103) 
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FIG. 11: The lowest energy quasiparticle state wavefunction density 
in the s-wave state, in the squar ^a)] circulai [(b)] rectangula j(c)| and 



elliptic (d) geometries. The density \u 3 -\ 2 + \v~j f 2 at each site j is rep- gulaj(c)|and elliptic (d) geometries. The density \uj 2 + 



FIG. 12: The lowest energy quasiparticle state wavefunction density 
in th e sp in-singlet p x + ip s state, in the squarc ^a)] circula ^b)] rectan- 

at each 



resented by the color of the dot at each lattice site. The highest den- 
sity is colored blue (black), and zero density is colored orange(light 
gray). It is clear in all cases that the lowest energy state is bound to 
vortex core at center of the geometry. Here the vortex core size was 
taken to be 0. With the nearest neighbor distance taken to be l/v3, 
the sq uarc^(a)| lattice patch has dimensions 12 x 12 (350 sites), the 
circulai [(b)| lattice patch has radius of 6 (262 sites), the rectangula (cj 
lattice patch has dimensions 10 x 15 (357 sites), and the elliptic (d) 
lattice patch has main axis of length 4 and 8 (222 sites). 



site j is represented by the color of the dot at each lattice site. The 
highest density is colored blue (black), and zero density is colored 
orange(light gray). It is clear in all cases that the lowest energy state 
is bound to vortex core at center of the geometry. Here the vortex 
core size was taken to be . With the nearest neighbor distance taken 
to be l/v3. the square (a) lattice patch has dimensions 12 x 12 (350 
sites), the circula^(b)]lattice patch has radius of 6 (262 sites), the rect- 
angul ai|(c)| lattice patch has dimensions 10 x 15 (357 sites), and the 
ellipticKdjlattice patch has main axis of length 4 and 8 (222 sites). 



where bj are bosonic operators and fj tCl are the fermionic op- 
erators used in our manuscript. As before i, j are used to de- 
note lattice sites, and the 2 greek letters a, (3 will be used to 
denote the spin indices f, |- The bosonic chemical potential 
is tuned to the value 3w - so that the bosonic band minimum 
is at zero energy, and po is the boson density per site. 

With no interaction between the boson and fermions, for 
small Ubb /w <C 1 the bosons will condense into a superfluid 
state. We assume we are deep in such a phase, and that the in- 
teractions with the fermions do not destroy the boson superflu- 
idity. We then use standard Bogoliubov theory to approximate 
the momentum space bose operators b M (q) « \/NoS(q) + 
6 M (q) where the second term is implicitly taken for only 
nonzero momentum. Here N is half the total number of 
bosons. The bosonic density operator then becomes p^(q) = 

£ k &t(k-q)& M (k) « at %) + ^[6t(-q) + ^(q)J- 

Using the Bogoliubov approximation we expand the bosonic 
terms of ( |103| l to quadratic order in the operators 6 M (q). Our 
goal is then to integrate out the bosonic degrees of freedom 



(the action is now Gaussian in the bosonic fields), and in this 
way find the effective Fermionic interactions that are induced. 

Taking only those terms in ( |103| l involving bosonic opera- 
tors, and performing a Fourier transformation we find 

n b =w ^(q) ( 3 <W - r A-(q)) Mq) 

+tfi//iVX;X> M (q)p /t (-q) 

q n 

+U bb /NY / J2 Mq) - Np Q S(q)} M-q) - Np d(q)} 

q |U 

(104) 

where we have introduced the fermionic density operator 
^(q) = Ek, Q fta ( k - q)/^^)- 11 is worth mentioning 
at this point that the boson density per site po = No /N where 
TV is the number of unit cells of the lattice. 

Assuming Nq is a macroscopic number, we can expand 
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Hb in powers of No. Keeping only the leading terms, we 
are left with a quadratic form in the bosonic operators. Next 
we apply a unitary transformation that diagonalizes the hop- 
ping term h (q) = ^e+ 10 / 2 [ai(q) + a 2 (q)] and 6 2 (q) = 

^e- l6j2 Mq) - a 2 (q)] where e+' e = 7(q)/|7(q)l (which 
implicitly depends on q). At this point it is worthwhile men- 
tioning that 9(— q) = — #(q), which is extremely useful in the 
detailed steps of our calculation that have been omitted here. 
After some rewriting of the Hamiltonian, we arrive at the re- 
markably separable form 



where for brevity we have written B = (q, r) . Integrating 
out the bosonic fields, we find 



q m 

-§Ei 



+'i > , K(q)M-q) + h.c] 



,+i0/2 



+ -L(Ft(q)e+«/ 2 -Ft(q)e- 



»e/2 



F 2 t(q) e -^ 2 ) 0l (q) 

) a 2 (q) + ^- c - 



+ ... 

(105) 



where we have introduced the coupling g — 2UbbPo, an d the 
band dispersions e\ = w (3 — |7(q)|)> £2 = w (3 + |7(q)|)- 
Note that all the momentum summations above formally ex- 
clude the q = mode. This will hold throughout the remain- 
der of this section, and so it will remain implicit. 

Next we employ a Bogoliubov transformation for each 
of the two bands (ai i2 ) separately. This is accom- 
plished by the transformation a M (q) = tt„(q).B(q) + 

/£ M (q)+e^(q)+g , , x _ 
2£ M (q) - 



Wu(q)-Bt(-q) with u„(q) 



The new operators £> (q) are canonical bosons, and the Hamil- 
tonian takes the form 

q fi 



+ (^ e +«/ 2 - ^V^ 2 ) (ua + «2) B 2 + /i.e. 



(106) 



where for the sake of brevity, we have omitted explicit men- 
tion of the q dependence of all the operators and functions 
such as v, u, E. 

In order to integrate out the bosonic fields, we must pass to 
a path integral formalism, taking into account the imaginary 
time derivative. Rewriting this term using the bosonic opera- 
tors (or complex fields) B^ (q) we find the action 
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where 



h \ = ^ (Fle+ ie /* + F 2 t e - ie / 2 ) ( Ul + Vl ) 



v^A 
v^A 



^ e +io/2_ F t e -ie/2y u2 + V2) 



(109) 



If we consider only low frequency effective interactions 
(uj n — » 0), then we can return to a Hamiltonian formulation of 
our problem with 



1 £ Ft(q) V(q)F,(q) . 



(HO) 



A 



The interaction vertex we have introduced V^(q) has com- 
ponents 



Vii(q) =V 22 (q) = U 2 f ^ 
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Note that Vij will be invariant under all the symmetries of the 
honeycomb lattice, including all lattice translations. 

Assuming the fermi level passes near the Dirac nodes 
of the honeycomb band structure (near half filling). The 
most significant low-energy excitations of the system then in- 
volve the fermionic operators with momenta near the Dirac 
nodes, at ±Q. The density fluctuation operator F^q) = 
J2k, a 4a( k - q)/ AlQ ( k ) wi H be significant only when k = 
±Q and k q = ±Q. This results in 3 possible regions for 
the exchange momentum: 

(i) q w 

(ii) q w 2Q = Q 

(iii) q w -2Q = Q . 

Expanding the function 7(q) around these 3 points, and as- 
suming the order of limits w 3> g S> wq 2 , we find 



K 
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N 2 
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Ag I2w ) 



f F 2 + F^F 1 



(112) 
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The case of non-interacting bosons g = 0, corresponds to 
a different limit than above w ^> wq 2 ^> g — 0. The Bo- 
goliubov spectrum becomes the same as the band structure 
E — > e, the Bogoliubov transformation parameters simplify 
to ii = 1 , v = 0, and the long-wavelength limit gives 



Vn(ci) = V 12 (ci) = U^\ 
6w q z 



(113) 



We find that in this limit the effective interactions are long- 
range in real space, and attractive. 

The boson interaction strength g controls the range for the 
effective interaction I ~ 1 j g. Despite us using a weak inter- 
action limit in the Bogoliubov theory, when assuming w ^> g 
we end up having the boson interaction coupling g dominating 
the nature of the effective interaction, with attractive q « 
interactions. The diagonal terms of the interaction vertex cor- 
respond to real space interactions between sites on the same 
sublattice. This type of interaction includes on-site interac- 
tions, which are expected to be the strongest of this kind. In 
contrast, the off-diagonal (the F^Fz term) terms of the inter- 
action vertex correspond to real space interactions between 
sites on different sublattices. The shortest range interactions 
in this class are nearest neighbor interactions. Since the off- 
diagonal and diagonal terms are comparable in magnitude, we 
expect an attractive effective nearest neighbor of comparable 
strength to that of the effective on-site interaction. 

From the features of the effective interaction vertex, we are 
lead to believe that the phenomenological model of Ref. |2T| 
may be an appropriate description of this system, which con- 
siders fermions on the honeycomb lattice, with only on-site 
and nearest neighbor interactions, parametrized by gi 2 re- 
spectively. For the range of parameters we find in the present 
work <72 < 0, and depending on the strength of Uff, the 
bare on site repulsion, we can have either positive or nega- 
tive sign of g\. Specifically, we can realize g\ < when 

4g ' 12w 



Uff « -^[is + mih and 5i > whent/^ » 



(± 

N 2 \4g 



1 

12w 



Uchoa et a/Pfind via a mean field analysis that the ground 
state may be a p + ip superconducting state for g 2 < and 
gi < and a mixed s-wave and p + ip superconducting state 
for g(2 < Oandgi > 0. Therefore it would seem that for strong 
bare fermion repulsion Uff, one should expect the p x + ip y 
phase. 



X. MAGNETIC FIELD SPLITTING 

Within the continuum Dirac theory we found zero modes 
in all the geometries we considered, topological protection oc- 



curs modulo symmetry mandated degeneracy (see section III I. 
The 4-fold degenerate zero modes we found in the two spin- 
singlet phases are protected as a result of the 4-fold degener- 
acy mandated by the SU(2) symmetries of the spin and the 
Dirac valley spinor. As first pointed out in ref. [19] a Zee- 
man field is found to split the spin-degenerate zero modes of 
the s-wave phase into Zeeman pairs, with a splitting propor- 
tional to magnetic field, at first order in perturbation theory. 



The magnetic field explicitly breaks the SU(2) spin symme- 
try, and therefore the zero modes can now split. In our ap- 
proach, it is easy to show this is an exact result, and that the 
zero mode states remain exact eigenstates of the system, albeit 
with a nonzero energy. 

Without a magnetic field the system is isotropic in the spin 
sector, and so we are free to choose the direction of the mag- 
netic field in spin space. It is convenient to choose the Zee- 
man field in the y-direction. The zero modes we found satisfy 
a v ip = a?p, and so we have precisely B<t v i/j = Bai\) splitting 
the zero modes. It is important to point out though, that while 
the mathematical spectrum of the BdG Hamiltonian has 4 zero 
modes splitting in to 2 positive and 2 negative, the physical 
excitations of the system consists only of the positive energy 
states (because of uj x Bct v uj x = Ba v , the BCS Hamiltonian 
identity ui x H.uj x = —H* still holds), and so there will be a 
doublet of lowest energy excitations, with E = B. 

Unlike a charged fermion superconductor, in a fermionic 
superfiuid there is no need for a magnetic field to create vor- 
tices. In any experimental cold atom apparatus, any magnetic 
field can be made extremely small. Therefore, it would be 
a great advantage to realize the spin-singlet phases we dis- 
cuss in this manuscript in a cold atom system, rather than a 
solid state system. If the magnetic field is too weak to destroy 
the fermion pairing, the only important effect of the mag- 
netic field is to introduce a Zeeman field coupling to the spin 
degrees of freedom, and the spectrum of vortex-core bound 
states can be manipulated. 

Using the Zeeman field splitting of the bound states, we 
now proceed to propose an experiment to probe whether zero 
modes exist in these systems (as found in the continuum Dirac 
theory) or not (expecting the zero modes to split slightly as 



found in the numerics of Section VIII ). The magnetic field 
allows us to control the low energy spectrum of the system, 
and this will possibly make it easier to identify in RF (low 
frequency) absorption measurements. If we assume all energy 
states of the experimental system are Kramers doublets, in the 
absence of zero modes, the lowest energy excitation (Ei) will 
be lowered when applying a magnetic field E = E\ — B 



(see Fig. 13 i. When zero modes exist, then when applying a 
magnetic field the lowest excitation energy will move up in 
energy (Eq — B) (see Fig. |T4]>. The lowest excitation energy 
in the system will then decrease with rising magnetic field in 
the absence of near-zero modes, but will increase if they exist 
in the system. This serves an experimental method to identify 
the existence of these states, which could easily be carried 
out in cold atom systems, but are perhaps more difficult in 
superconducting solid states systems. 



XI. DISCUSSION AND CONCLUSION 

In this manuscript we have explored whether topological 
zero modes exist in a number of possible fermionic conden- 
sate phases on the honeycomb lattice. We examined 2 spin- 
singlet phases both of which are fully gapped in the entire 
Brillouin zone. We have found that 4-fold degenerate topo- 
logical zero modes exist within the continuum Dirac theory, 
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FIG. 13: Influence of the Zeeman field splitting on the excitation 
spectrum in the case that vortex core zero modes are absent 



E=E +B 



E=E, 



first excited mode 



E=E -B 



zero modes 



E=0 



E=B 



FIG. 14: Influence of the Zeeman field splitting on the excitation 
spectrum in the case that vortex core zero modes exist 



for these two phases. We have done this by explicitly solving 
for zero modes bound to vortex cores, at sample edges, and 
in SNS junction geometries. In all cases, the edge state and 
vortex core calculations agree completely, with the same de- 
generacy of zero modes being found, but in the SNS junction 
geometry an extra accidental symmetry doubles the number 
of zero modes from 4 to 8. 

With an even degeneracy of vortex core zero modes, the 
majorana zero modes are not compelled to pair into fermionic 
degrees of freedom between spatially remote vortices, but 
rather locally at each vortex core. The natural mechanism 
for entanglement of vortex pairs is therefore lost, and no non- 
Abelian statistics between vortices should appear. 

The topological zero modes existence crucially depends on 
the emergent low energy SU (2) symmetry of the Dirac valley 
spinor structure in the Dirac theory. This symmetry does not 
strictly hold in the original lattice model, and this brings up the 
possibility that the seemingly protected zero modes found in 
the Dirac theory, are split in the full honeycomb lattice model. 
While corrections to the Dirac theory give an unclear picture 
of the fate of the zero modes when the effective low energy 
symmetry is broken, the numerical diagonalization we have 
performed on the lattice model confirm that indeed this is the 
case - the zero modes are split. In the present context simply 
using the continuum Dirac model is inaccurate - even when 
the vortex structure is slowly varying on the scale of the lat- 
tice, the Dirac theory is still only approximate, and in fact the 
zero modes are not topologically protected. 

We have also discussed the Zeeman field splitting of the 
vortex core bound states in this phases. We suggested an 
experiment taking advantage of this splitting to ascertain 
whether zero modes exist or not in this system, by tracking 



how the excitation energies change when modifying the Zee- 
man field. 

In order to realize the experiment we propose, one first 
needs to create a condensate on the honeycomb lattice. It may 
be possible to realize a superconducting state in graphene by 
superconducting leads inducing electron pairing via the prox- 
imity effect. Another possibility we have discussed, is form- 
ing fermion condensates in cold atom Bose-Fermi mixtures. 
The latter however will most probably require the fermions 
to be cooled down to very low temperatures compared with 
the energy scales of optical lattices, which is challenging in 
current experiments. A cold atom gas is perhaps the ideal re- 
alization of the condensate for our proposed experiment, since 
magnetic fields are not involved in the forming of vortices in 
the first place, and so can be freely manipulated, without af- 
fecting the condensate or the vortices too much. 

In the context of the zero modes in the "ordinary" p x + ip y 
state, it has already been suggested to probe the bound state 
spectrum by RF absorption 37 , and STM measurements 3 -^. The 
same tools could be used to probe the bound state spectrum 
in the phases we discuss here. As opposed to Refs. 1371381 in 
the experiment we propose one would be looking for how the 
spectrum moves when changing the magnetic field, rather than 
simply looking at a static spectrum. It is sometimes more easy 
to notice something that is moving, rather than stationary, and 
so it may prove easier to detect the spectrum shifts. 

The reaction to magnetic field of the bound states is not lim- 
ited to vortex cores - it could be discernible in edge states as 
well if the sample is small enough that the discreteness of the 
energy levels bound to the edge becomes evident. The edge 
state spectrum in all the phases we consider here is always 
linear in the transverse momentum E ~ qv, with some ef- 
fective velocity v. For a finite system, the momentum will be 
quantized, q = {I + 7), and the bound states have a level 
spacing of ~ ■ If tne Zeeman splitting as well as the ther- 
mal energy scale are smaller than this level spacing, the effect 
we describe here is in principle observable. In practice, the 
experimental probe must be sensitive enough to probe these 
small energy scales. 

Finally, our present work has forced us to generalize some 
ideas that were understood and developed in the context of the 
"ordinary" p x + ip y state. Topological zero modes were pre- 
viously understood to be topologically protected only if single 
zero modes existed (in a unit vorticity vortex). We have gen- 
eralized this view of topological protection to accommodate 
symmetry mandated degeneracy of quasiparticle excitations, 
for which the phases of the continuum Dirac theory we con- 
sidered are examples . Our conclusion is that zero modes can 
be protected to all perturbations that preserve the symmetries, 
and as such are protected by the combination of symmetry 
and topology. If the symmetries are explicitly broken by a 
perturbation, then the zero modes may split. This is precisely 
the reason the Dirac theory and the precise honeycomb lat- 
tice model differ, and also the reason for the Zeeman splitting 
of the zero modes in the Dirac theory. This synergistic protec- 
tion, while clearly more fragile than the topological protection 
of a single zero mode may prove important in understanding 
many other physical systems beyond those we discuss here. 
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We have also generalized the known connection between 
topological zero modes bound to vortex cores and at sample 
edges in the well studied "ordinary" p x + ip y state. We have 
shown that in quite general settings, with the possibility of 
symmetry mandated degeneracy included, zero modes bound 
to vortices and edges should be identified. 

Our analysis was limited to a number of presumed pair- 
ing states for fermions on the honeycomb lattice, but similar 
phenomena may be uncovered in other states involving the 
honeycomb lattice. In particular Refs. 39 40 have discussed 
non-condensate models on the honeycomb lattice with vor- 
tices possessing zero mode bound states in their core. In both 
cases only a vortex calculation was carried out, and given our 
general observation that vortex core bound zero modes should 
be identified with edge state zero modes, we expect these zero 
modes to appear at sample edges in the models of Refs. [391401 
In Ref. 40 the authors find there is a single zero mode bound to 
the vortex core. Another related model which exhibits topo- 
logical zero modes at sample edges is the Kane-Mele model 41 , 
which both in the precise lattice model, as well as in a contin- 
uum limit 42 , exhibits edge state zero modes (in the continuum 
case, for an armchair boundary - the zigzag boundary suffers 
from the same problems we pointed out in sectionfVl}. Finally, 
we mention a very recent publication 4 ^ finding zero modes 
bound to vortices in a bilayer-graphene exciton condensate. 
As in our case, the zero modes turn out to split in the precise 
lattice model. We suspect that other interesting possible states 
of matter on the honeycomb lattice geometry exist, as well as 
in 3-dimensional geometries that supply the common ingre- 
dient in all these models - the Dirac nodes in the lattice band 
structure. The physics of a 3-dimensional version of the Kane- 
Mele model 44 is realized in Bii-^Sb^, as recently probed in 
Ref. |45] and following this work Ref. |46]has suggested that 
majorana fermion zero modes should appear at the interface 
between a topological insulator and an s-wave superconduc- 
tor. 
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APPENDIX A: SPINLESS P+IP CONDENSATE 

In this appendix we analyze the spinless p x + ip y phase 
mentioned briefly in the main text. We will find that this phase 
in some geometries will posses zero modes, but these are bulk 
states rather than bound states. 



1. SNS junction 

In this subsection we will analyze the SNS junctions in the 
spinless p x + ip y phase. We consider here only wavefunc- 



tions that are uniform in the direction parallel to the SNS junc- 
tion boundaries, and find in stark contrast to the spin-singlet 
phases, that no zero modes exist. The steps of our analysis fol- 
low closely those of the SNS calculations for the spin-singlet 
phases, and so we will describe our calculations in minimal 
detail. 

With a (piecewise) uniform pairing function, combining the 
kinetic energy pi) and the spinless p x + ip y pairing term p4| 
yields a BdG equation of the form 



7-tBdG^ = Etp 
\iio z — ivD + 



c Ar] z D [u x cos(<p) + uj v sin(</>)] ^ 



(Al) 



where we have dropped the 2 subscript from both the order 
parameter phase and magnitude, to avoid clutter. 

We consider only states that are uniform in the direction 
parallel to the SNS junction boundaries, and use the same uni- 
tary transformation U(a) to rotate the angle between the SNS 
boundaries and the y-axis a — ► 0, in the operator D(a) ap- 
pearing in both the kinetic energy and the pairing term. All 
other parts of the BdG Hamiltonian remain invariant, and as 
long as we ignore the quadratic correction ( |25) , we can simply 
set a — 0. The BdG equations reduce to 

fiuj z -E-ivri x d x +T x Aid x [u> x cos(</>) + uj y sm((p)] tp = . 

(A2) 

As before, we will work in the London gauge which we can 
get by applying the unitary transformation O = er 1 ^^'. We 
are left with 



(iu z — E — ivifd x + r x Aid x 



4> = o 



(A3) 



from which it is clear that we can choose solutions that are 
eigenstates of both r\ x and t x , such that r) x ip = rjtp and t x tp = 
Ttp. The BdG equations then can be reorganized in the form 
id x tp — Atp with 



1 / vr)(E - n) A(E + h)t 

-A 2 V A(E-h)t vr](E + [i) 



(A4) 



The eigenvalues of the matrix 

V /A 2 E 2 + n 2 (v 2 



A are vErj ± 



A 2 ), and since we expect v 3> A, 
the eigenvalues will in general be real numbers. As a conse- 
quence, we can only have solutions of tp that are exponentials 
of purely imaginary arguments. As a result, no bound states 
can appear - these require some exponential decay of the 
wavefunction. 



2. Edge geometry 

We will now address the edge geometry in the spinless p x + 
ip y phase. Starting form the BdG equation \A\\ , with the 
phase chosen as <p = 0, and assuming a transverse momentum 
q, such that tp = e lqv ip(x). The reduced BdG equations are 



[1UJ 



ivD + T x Ar/ z Duj x ri y - E 



tp = 



(A5) 
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with D = rfd x + u> z rfd y . Reorganizing the BdG equations 
yields 

\[iu z -E + vqr/ v uj z - qAr x rj z uj v } ip = i [vrf - At x uj x ] d x ip 

(A6) 

Using the identity [vrj x — At x oj x ] = 
V 2^ A 2 [vrf + At x uj x ], we bring the equation to the 
form d x ip — Aip, with 



1 



(vr] x + At x lo x ) {[iu z - E) + iquj z rj z 



A 



(A7) 

We are free at this stage to choose solutions that are eigen- 
states of T x ip — Tip, so we simply replace t x — ► r. Fur- 
thermore, it is useful at this point to denote z = ^ff^a » 

e = — and k = *?V" 2 ~ A2 a u f which are dimensionless 
quantities. We can also assume without loss of generality that 
H > 0, v > A, and E > 0. The BdG equations now become 
d z ip — A with 



.4 



1 



{yrf + At x u x ) (uj z - e) + iku z rj z 



Vv 2 - A 2 

The 4 eigenvalues of the matrix iA in this notation are 

±A„ = 

..2 L A 2 



± 



1 - fc^ + 



J(v 2 - A 2 ) + A 2 e 2 

^ 



(A8) 



1 1/2 



(A9) 



where r\ = ±1. 

We are interested in exploring zero modes, so at this point 
we set e = 0, to identify which eigenvalues can give solu- 
tions that are exponentially decaying in the z > region. The 
eigenvalues of A become ±A, ; = ±i\/l — k 2 . For exponen- 
tially decaying zero modes, we must therefore have \k\ > 1. 
In this regime, q > ^ V ^_ A2 ■ It i s noteworthy that the nodes 
in the bulk spectrum for this phase occur in the continuum 

Therefore, these 



theory at precisely q = (-s/3, l) 
bound states may be identified with the bulk zero modes. 



V / 1 ,2_ A 2 ' 



3. Vortex geometry 

In this subsection we turn to explore whether zero modes 
exist bound to vortex cores in the spinless p x + ip y phase. 

We start with the BdG equation [Ho + H3] ip = Eip with 
Tio, 3 from pTj ), and ( |24] i. We choose the eigenstates to satisfy 
r x ip = rip, as in the SNS and edge geometries. We will model 
the vortex by assuming the form iA2 = A(r)e +l ^, with A(r) 
real (as for the spin singlet p x + ip y case, the convention is 
different from earlier sections so that we can use the same 
ansatz for the polar angle dependence as for the s-wave case). 
Also, since it will prove convenient, we will assume that the 
order parameter radial profile is piecewise uniform. 

As in previous subsections, we find the (p dependence can 
be eliminated from the zero-mode problem by choosing the 



wavefunction form 

iP(r,(P) = e 4 ^ (u 1 (r) 1 e^u 2 (r),v 1 (r),e-^v 2 (r)) T . 

(A10) 

The reduced ODEs then involve only the radial coordinate, 
and can be cast in the form d r ip = Aip where 



A 



1 



1 



lw Z T) Z 



r- 2 
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i] z + ivi] x (E - i±uj z ) 



v 2 - A 2 

V —iu x T) v + AT(Eicu x - nu v ) 
2r 



(All) 



For the purpose of showing that no bound state zero modes 
exist, it will suffice to consider the asymptotic limit (r — > 00) 
alone. We neglect all terms in A that have a factor 1/r 



1 



A = -^-r^ [iv V x (E - ^u z ) + A T {Eiu x - ^ v )} . 
v — A z 

(A12) 

Setting E — 0, we find A has the eigenvalues ±i^==j, 
which are purely imaginary. Therefore no bound states with 
zero energy are allowed! Thus, we conclude that no zero 
modes exist at the vortex core in this phase. 

This particular calculation shows us that having a Dirac 
equation structure in the BdG equations, is not a sufficient 
condition for topological zero modes to be present. 



4. Numerics 

In the previous sections of this appendix, we found that 
the edge state geometry can support some zero modes with 
a wavefunction concentrated at the at the edge. In contrast, in 
the vortex calculation found no bound states zero modes. As 
we argued in this paper, in a fully gapped system, we expect a 
general correspondence between the edge state spectrum and 
the vortex core bound state spectrum. This expectation does 
not hold here, presumably due to the fact that this phase is not 
fully gapped. 

To verify that the vortex calculation result is correct (it uses 
the approximate continuum description) we employed the nu- 



merical methods of section VIII Using the precise lattice pair- 
ing function for the spinless p x + ip y phase, with a square lat- 
tice patch of 1824 lattice sites, we calculated the low energy 
spectrum for the vortex state, with two different vortex core 
sizes R — 0, v5. We set the parameters | A| = 0.5 and \x = 
0.4. The lowest energies divided by the de Gennes scale are 
= ±0.0270939, ±0.0558193, ±0.0765734 . . ., for R = 



E 

A 2 



and 



±0.0234966, ±0.0350545, ±0.0630661 ... for 



R = \/5. The results in both cases are similar - we find low 
energy states exist, far below the de Gennes energy scale, but 
all the states with energy below the de Gennes scale, are delo- 
calized bulk states and not concentrated near the vortex core. 
This result would indicate that the vortex calculation and edge 
state calculations are not at odds - the zero modes found in the 
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edge state calculation are related to the bulk low energy states that appear due to the nodes in the pairing function. 
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